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Physical Processes in Protoplanetary Disks 

Philip J. Armitage 


Abstract This review, based on lectures given at the 45th Saas-Fee Advanced 
Course “From Protoplanetary Disks to Planet Formation”, introduces physical pro¬ 
cesses in protoplanetary disks relevant to accretion and the initial stages of planet 
formation. After reprising the elementary theory of disk structure and evolution, 
I discuss the gas-phase physics of angular momentum transport through turbulence 
and disk winds, and how this may be related to episodic accretion observed in Young 
Stellar Objects. Turning to solids, I review the evolution of single particles under 
aerodynamic forces, and describe the conditions necessary for the development of 
collective gas-particle instabilities. Observations show that disks are not always ra¬ 
dially smooth axisymmetric structures, and I discuss how gas and particle processes 
can interact to form observable large-scale structure (at ice lines, vortices and in 
zonal flows). I conclude with disk dispersal. 


1 Preamble 

The objective of this review is to introduce the physical processes in protoplanetary 
disks that are relevant to protostellar accretion and the initial stages of planet forma¬ 
tion. Protoplanetary disks, as well as being interesting objects of study in their own 
right, are also simultaneously an outcome of star formation and initial conditions 
for planet formation. As such, we need to understand the evolution of both the dom¬ 
inant gaseous component and the trace of solid material that is critical for planet 
formation. Much interesting complexity occurs due to interactions between the two. 

The review is organized around three motivating questions; how do protoplane¬ 
tary gas disks evolve over time, how are solids transported and concentrated within 
the gas, and how do gas-phase and solid processes interact to form structure within 
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disks? Protoplanetary disks are long-lived entities, and I begin in plby outlining 
their equilibrium structure. Disk evolution is described in fj3] and f|4| first follow¬ 
ing the classical approach in which the origin of angular momentum transport is 
unspecified, and subsequently in a more modern presentation where angular mo¬ 
mentum transport and loss processes are ascribed to specific fluid instabilities. lj5] 
discusses candidate theoretical explanations, some of them directly related to angu¬ 
lar momentum transport processes, for episodic accretion outbursts. In <|6]the focus 
switches to solids, and I review how single solid particles settle, drift, diffuse and 
concentrate relative to the gas. fjTJthen describes how these processes can give rise 
to structure within the disk on various scales, either at transient “particle traps” or 
at persistent locations such as ice lines where the disk structure varies rapidly. Fi¬ 
nally, previews what is known about the processes, foremost amongst which is 
photoevaporation, that can disperse protoplanetary disks. 

The lectures, even distended as in this review, could not touch upon all of the 
physics that an aspiring researcher in the field might require. In an effort to cover as 
much as possible — and to accommodate the diverse backgrounds of participants 
(and readers) — the material is presented with varying degrees of detail and rigor. 
For most topics I begin with a self-contained discussion of essential material that 
would often be assumed as background in papers or talks. I then discuss the under¬ 
lying physics of more recent results, and give entry points to the relevant literature. 
The gray boxes in the text identify open issues that interest me. 

One caution is in order. Quantities are generally defined and labeled following 
conventions in the literature, to make it easier for the reader who needs to fill in 
missing details or to explore further. Given the broad scope of the review there is 
considerable overloading of notation, with the same symbols being used to represent 
unrelated quantities in different sections. Take care! 


2 Disk structure 

The lifetime of protoplanetary disks is observed to be of the order of several Myr 
1 13811147 1. which equates to millions of dynamical times in the inner disk and thou¬ 
sands of dynamical times in the outer disk at 100 AU. To a first approximation 
we can treat the disk as evolving slowly through a sequence of axisymmetric static 
structures as mass accretes on to the star and is lost through disk winds, and our first 
task is to discuss the physics that determines those structures. Quantities that we are 
interested in include the density p(r,z), the gas and dust temperatures T(r,z ) and 
Td(r,z), the chemical composition, and the ionization fraction. The importance of 
these for observations and for planet formation is self-explanatory, with the possible 
exception of the ionization fraction which determines how strongly magnetic fields 
couple to the gas. The density of solid particles (“dust”) p f / is also important, but we 
will defer saying much about that until we have discussed turbulence, radial drift, 
and the aerodynamic coupling of solids and gas. 
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Fig. 1 The geometry for calculating the vertical hydrostatic equilibrium of a non-self-gravitating 
protoplanetary disk. The balancing forces are the vertical component of stellar gravity and the 
vertical pressure gradient. 


2.1 Vertical and radial structure 


2.1.1 Vertical structure 


The vertical profile of gas density in protoplanetary disks is determined by the con¬ 
dition of hydrostatic equilibrium. The simplest case to consider is an optically thick 
disk that is heated by stellar irradiation, has negligible mass compared to the mass 
of the star, and is supported by gas pressure. We can then approximate the optically 
thick interior of the disk as isothermal, with constant sound speed c s and pressure 
P = pc 2 . The sound speed is related to the temperature via c 2 = k^T / pmn, where ks 
is Boltzmann’s constant, m/y is the mass of a hydrogen atom, and where under nor¬ 
mal disk conditions the mean molecular weight jl ~ 2.3. In cylindrical co-ordinates, 
the condition for vertical hydrostatic equilibrium (Figure [TJ is, 

dP GM * ■ a 

-^ = -pg z = — f 2-^sm0p, ( 1 ) 

where M* is the stellar mass. For z<r, 

GM* , 

fc= (r 2 +z 2 )3 /2 ^^ (2) 


where f2 K = GM*/r 3 is the Keplerian orbital velocity (here defined at the mid¬ 
plane, later in 12.1.2 we will need to distinguish between the mid-plane and other 
locations). Equation |T|i then becomes. 


2 ([P 

C * d" 




which integrates to give. 


P(z) =Poexp [~z 2 /2h 2 ] , 


(3) 

(4) 


where po is the mid-plane density and we have defined the vertical scale height 
h = c s /£ 2k- Because the effective gravity increases with height (and vanishes at the 
mid-plane) this standard disk profile is gaussian, rather than exponential as in a thin 
isothermal planetary atmosphere. A consequence is that the scale over which the 
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Fig. 2 Simple models for the vertical density profile of an isothermal disk, in units of the disk 
scale height h. The solid blue line shows the gaussian density profile valid for the dashed 

blue line shows the exact solution relaxing this assumption (the two are essentially identical for 
this disk, with h/r = 0.05). The red dashed curve shows a fit to numerical simulations that include 
a magnetic pressure component G50). 


density drops by a factor of e gets smaller with z; loosely speaking disks become 
more “two dimensional” away from the mid-plane. Defining the surface density 
E = / pdz, the central density is, 


1 E 


(5) 


Up to straightforward variations due to differing conventions (e.g. some authors de¬ 
fine h = v2Ci/X2K) these formulae define the vertical structure of the most basic 
disk model (isothermal, with a gaussian density profile). For many purposes, espe¬ 
cially if one is mostly worried about conditions within a few h of the mid-plane, it 
may be a perfectly adequate description. 

The most obvious cause of gross departures from a gaussian density profile is a 
non-isothermal temperature profile. If the disk is accreting, gravitational potential 
energy that is thermalized in the optically thick interior will need a vertical tem¬ 
perature gradient AT/Az < 0 in order to be transported to the disk photosphere and 
radiated. We will return to this effect in ' \2.2\ and Sj3] after assessing the other as¬ 
sumptions inherent in equation (j^J. We first note that the simplification to z r is 
convenient but not necessary, and that we can integrate equation G without this 
assumption to give, 




Physical Processes in Protoplanetary Disks 


5 


p = po exp 


h 2 


(l+r/A 



( 6 ) 


Protoplanetary disks are geometrically thin, however, with h/r ss 0.05 being fairly 
typical. In this regime, as is shown in Figure [2] departures from a gaussian are 
negligibly small. We only really need to worry about the full expression for vertical 
gravity when considering disk winds, which flow beyond z ~ r. 

What about the contribution of the disk itself to the vertical component of grav¬ 
ity? Approximating the disk as an infinite sheet with (constant) surface density E, 
Gauss’ theorem tells us that the gravitational acceleration above the sheet is inde¬ 
pendent of height, 

g z =2nGE. (7) 

Comparing this acceleration with the vertical component of the star’s gravity at 
Z = h, we find that the disk dominates if, 


E 


M h 
2 k r 3 


( 8 ) 


Very roughly we can write the disk mass Mdisk ~ Kr 2 E, which allows us to write the 
condition for the disk’s own gravity to matter as, 


4^disk 1 / h \ 
M* > 2\r) 


(9) 


For h/r = 0.05, a disk mass of a few percent of the stellar mass makes a non- 
negligible change to the vertical structure, and such masses are not unreasonably 
large. As we will see in fjd] however, when disk masses /Wdisk/AC ^ } x / r are encoun¬ 
tered we tend to have bigger fish to fry, as this is also the approximate condition for 
the onset of disk self-gravity, the formation of spiral arms, and substantial departures 
from axisymmetry. 

Magnetic pressure Pb = B 2 /8k is likely to impact the vertical density profile, 
at least for z h (here and subsequently, we use units such that B is measured in 
Gauss). No simple principle for predicting the strength or vertical variation of the 
magnetic field is known, so we turn to numerical simulations for guidance. Hirose 
& Turner Q30) completed radiation magnetohydrodynamic (MHD) simulations of 
the protoplanetary disk at 1 AU, taking the surface density E = 10 3 g cm~ 2 and 
adopting as stellar parameters A7, : =0.5 M 0 , r e ff = 4000 K, and IP = 2 R .. Their 
simulations included Ohmic diffusion but ignored both ambipolar diffusion and the 
Hall effect (see f[4] for further discussion of these processes). They found that the 
model disk was gas pressure dominated near the mid-plane, but that the atmosphere 
(or corona) was magnetically dominated. An empirical fit to their density profile is 

8333 , 

P = [exp (-r/2/i 2 ) +eexp(-\z\/kh)\ , (10) 

with e ~ 1.25 x 10~ 2 and k ~ 1.5. This fit is shown in Figure]^ Magnetic pressure 
beats out gas pressure for kl > Ah, leading to a low density exponential atmosphere 
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that is much more extended than a standard isothermal disk. Even if the atmosphere 
itself gives way to a disk wind at still higher altitudes (as suggested by other simu¬ 
lations), these results suggest that observational probes of conditions near the disk 
surface may be sampling regions where the magnetic field dominates. 


The impact of magnetic pressure on the vertical structure of protoplanetary 
disks may be larger than the preceding discussion implies. Simulations of the 
inner disk (r ~ 1 AU) that include the Hall effect lj 1951 result in the genera¬ 
tion of strong azimuthal magnetic fields whose magnetic pressure may exceed 
that of the gas within one scale height, or even at the mid-plane. This implies 
that there could be regions of the disk where the conventional estimate of the 
scale height based on the temperature is significantly in error. Observational 
constraints on such models would be valuable, as would a better analytic un¬ 
derstanding of what determines the vertical profile of B in disks. 


2.1.2 Radial structure 

The radial run of the surface density cannot be predicted from considerations of 
static disk structure; it is either an observational question 0 or one to address using 
time-dependent models (fj3j^| Instead, we consider the azimuthal velocity v l? . If the 
disk is static (and even if it is slowly evolving) the azimuthal component of the 
momentum equation. 


<9v 1 

+ (v- V)v =-VP- V<£> 

dt v ' p 


can be written in the mid-plane as, 

_ GMh, | 1 dP 
r r 2 p dr 


(ID 


( 12 ) 


Here P is the pressure and all quantities are mid-plane values. Let’s start with an 
explicit example of the consequences of this force balance in protoplanetary disks. 
Consider a disk with Zap 1 and central temperature T c r l2 . We then have 
c s r -1 / 4 , p r~ 9 / 4 and P r -11 / 4 . Substituting into equation (12 1 yields. 


= VK 1 



(13) 


1 The minimum mass Solar Nebula (MMSN) {3471ll42t . an approximate lower bound for the 

amount of disk gas needed to form the planets in the Solar System, can be useful as a reference 

model despite its tenuous connection to actual conditions in the disk at the time of planet formation. 

The MMSN has a gas surface density profile E(r) = 1.7 x 10 3 (r/A[/)~L 2 g cm~ 2 . 
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From this we deduce, 

• The deviation from strict Keplerian rotation, vk = \JGM, /r is of the order of 

{h/r) 2 . 

• Its magnitude is small. For a disk with h/r = 0.03 at 1 AU, the difference be¬ 
tween the disk azimuthal velocity and the Keplerian value is about 0.25%, or, in 
absolute terms |v^ — vk| — 70 m s . 

When we come to discuss the evolution of particles within disks (Sj6|, it will turn 
out that this seemingly small effect is of paramount importance. Particles do not 
experience the radial pressure gradient that is the cause of the mismatch in speeds, 
and as a result develop a differential velocity with respect to the gas that leads to 
aerodynamic drag and (usually) inspiral. Because this process is so important it is 
worth studying not just the magnitude of the effect but also its vertical dependence. 
To do so, we follow Takeuchi & Lin 02211 and consider an axisymmetric vertically 
isothermal disk supported against gravity by gas pressure. The vertical density pro¬ 
file is then gaussian (equation [4]> in equilibrium (equation [TT| we have. 


rQ 2 = 


GM* 


1 dP 


(rl+z 2 ) 312 ^ P dr ' 


(14) 


We distinguish between the gas angular velocity, X2 g (r,z), the Keplerian angular 
velocity £2^{r,z) = GM,/(r +z 2 ) 3 / 2 , and its mid-plane value ^K.mid- The disk is 
fully specified by the local power-law profiles of surface density and temperature. 


E oc r ~ y (15) 

Tc^r-P, (16) 


with 7=1 and J3 = 1/2 being typically assumed values. Evaluating dP/dr using 
equation Q with h = h(r) allows us to determine the equilibrium gas angular ve¬ 
locity in terms of the mid-plane Keplerian value. 


i2g — TlK.mui 


1 - U£) ( p+2y+3+ ^h 


(17) 


Provided that the temperature is a locally decreasing function of radius (ji > 0), the 
sense of the vertical shear is that the gas rotates slower at higher z. Like the sub- 
Keplerian mid-plane velocities, the magnitude of the shear is only of the order of 
{h/r) 2 , but this small effect may nevertheless be detectable with ALMA data B297L 
For particle dynamics, what matters is the difference between the gas velocity and 
the local Keplerian speed. To order z 2 /r 2 , the vertical dependence of the Keplerian 
velocity is, 

^K^ K , mid (l~^. (18) 

This function decreases faster with height than f2 g , so the difference between them, 
plotted in Figure [3] switches sign for sufficiently large z, 







Philip J. Armitage 


\ 

Sh 


•D 

£ 

Sd 

a 

\ 

a 

<3 



Fig. 3 The variation in angular velocity with height in the disk. In blue, the angular velocity of gas 
relative to the mid-plane Keplerian value, (fl g — f2K,mid)/^K.mid x ( r /h)~- In red, the difference 
between the Keplerian angular velocity and the local gas angular velocity, (i2 g — /flionid x 
(r/h) 2 . The assumed disk has £ °c and T c 2 (solid curves) or radially constant T c (dashed 

curves). 


£2 g — 


1 ( h 

4 \ r 


j3+2y+3 + (/3-3)^ 


n 


K.mid- 


(19) 


Particles that orbit the star at the local Keplerian speed move slower than the gas 
near the mid-plane (and thus experience a “headwind”), but faster at high altitude. 
For typical parameters, the changeover occurs at about z « 1.5/;. 

The orbital velocity will also deviate from the point mass Keplerian form if the 
disk mass is sufficiently high. The gravitational potential of a disk is not that of 
a point mass (and does not have a simple form for realistic disk surface density 
profiles), but for an approximation we assume that it is. Then the modified Keplerian 
velocity depends only upon the enclosed disk mass. 


V K - ' ; K 



1/2 


( 20 ) 


For disk masses M^sk ~ 10 2 A 7 /, the effect on the rotation curve is of comparable 
magnitude (but opposite sign) to the effect of the radial pressure gradient. 
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Fig. 4 The setup for calculating the radial temperature distribution of an optically thick, razor-thin 
disk. We consider a ray passing that makes an angle 6 to the line joining the area element to the 
center of the star. Different rays with the same 6 are labeled with the azimuthal angle <j >; ip = 0 
corresponds to the “twelve o’clock” position on the stellar surface. 


The simplicity of the calculation of the deviations from Keplerian rotation 
due to pressure gradients in the disk depends on the disk being axisymmet- 
ric. How axisymmetric (and planar) are real disks? In most cases we don’t 
know observationally, and it’s obvious that a disk formed from the collapse of 
a turbulent molecular cloud core might well start out significantly eccentric 
and warped. The rate of decay of eccentric or warp perturbations remains a 
subject of active research 1 45 1. It’s thus worth remembering, especially when 
interpreting precise kinematic observations, that un-modeled warps or eccen¬ 
tricities as small as e ~ ( h/r ) 2 could be significant. 


2.2 Thermal physics 

We seek to determine the temperature of gas and dust as Our first task is 

to calculate the interior temperature of a disk heated solely by starlight. This is 
straightforward. At most radii of interest the dust opacity is high enough for the 
disk to be optically thick to both stellar radiation and to its own re-emitted radiation, 
which hence has a thermal spectrum. It is then a geometric problem to work out how 
much stellar radiation annulus of the disk intercepts, and what equilibrium T results. 
We will then consider the temperatures of gas and dust in the surface layers of disks. 
These problems are trickier. The surface layers are both optically thin and of low 
density, so we have to account explicitly for the heating and cooling processes and 
allow for the possibility that the dust and gas are too weakly coupled to maintain the 
same temperature. We defer until f|3]the question of how accretion heating modifies 
these solutions. 

A disk whose temperature is set by stellar irradiation is described as “passive”. 
The model problem is a flat razor-thin disk that absorbs all incoming stellar radiation 
and re-emits it locally as a blackbody. We seek the temperature of the blackbody disk 
emission as f(r). Modeling the star as a sphere of radius R ,., and constant brightness 
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/*, we define spherical polar coordinates such that the axis of the coordinate system 
points to the center of the star (figure |4|. The stellar flux passing through a surface 
at distance r is, 

F = J /* sin 9 cos (j)d £2, (21) 

where dQ represents the element of solid angle. We count the flux coming from the 
top half of the star only (and later equate that to radiation from only the top surface 
of the disk), so the integral has limits, 

- n/2 < <j> < 7t/2 

0 < 0 < sin -1 • (22) 

Substituting dQ = sin ddQd<f>, the integral evaluates to. 



A star with effective temperature 7) has brightness /* = (1/k)gT*, with a the 
Stefan-Boltzmann constant. Equating F to the one-sided disk emission o7~[j sk the 
temperature profile is. 



The exact result is unnecessarily complicated. To simplify, we expand the right hand 
side in a Taylor series for (R*/r) <C 1 (i.e. far from the stellar surface) to obtain, 

Tdisk - (25) 

as the power-law temperature profile of a thin, flat, passive disk. This implies a 
sound speed profile, c s r 3/8 , and a disk thickness (h/r) °c r |/x . We therefore 
predict that the disk becomes geometrically thicker (“flares”) at larger radii. We can 
also integrate equation ( f24| > exactly over r, with the result that the luminosity from 
both side of the disk sums to L^k/L* = 1 /4. 

A more detailed calculation of the dust emission from passive disks requires 
consideration of two additional physical effects cm First, as we just noted, the 
disk thickness as measured by the gas scale height flares to larger radii. If dust is 
well-mixed with the gas — which may or may not be a reasonable assumption — 
a flared disk intercepts more stellar radiation in its outer regions than a flat one, 
which will tend to make it flare even more strongly. We therefore need to solve 
for the self-consistent shape of the disk that simultaneously satisfies hydrostatic 
and thermal equilibrium at every radius. This is conceptually easy, and the slightly 
messy geometry required to generalize the flat disk calculation is clearly described 
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by Kenyon & Hartmann 11 171 II . Second, small dust grains that are directly exposed 
to stellar irradiation (i.e. those where the optical depth to stellar radiation along 
a line toward the star T < 1) emit as a dilute blackbody with a temperature higher 
than if they were true blackbody emitters (EZD . The reason for this is that small dust 
grains, of radius s, have an emissivity e = 1 only for wavelengths A < 2ns. At longer 
wavelengths, their emissivity declines. The details depend upon the composition 
and structure of the dust grains, but roughly the emissivity (and opacity k) scale 
inversely with the wavelength. In terms of temperature. 


£ = 



(26) 


with /3 = 1. A dust particle exposed to the stellar radiation held is then in radiative 
equilibrium at temperature T s when absorption and emission are in balance. 


4 nr 2 


ns~ 


oT s 4 e(T s )4ns 2 . 


(27) 


The resulting temperature. 




1 

jrj4 



(28) 


exceeds the expected blackbody temperature by a substantial factor if e 1. 

An illustrative analytic model that incorporates these effects was developed by 
Chiang & Goldreich ||75| . They considered a disk with a surface density profile E = 
10 3 (r/l AU)~ 3 / 2 g cm" 2 around a star with M* = 0.5 Af 0 , T* = 4000 K and = 
2.5 Rq. Within about 100 AU, their solution has half of the bolometric luminosity 
of the disk emitted as a blackbody at the interior temperature, 

i '‘ s = 150 (txc)" ,, k ' (29) 

with equal luminosity at each radius emerging from a hot surface dust layer at. 



The Chiang & Goldreich solution is a two-layer approximation to dust continuum 
radiative transfer for a passive, hydrostatic disk. Approximations in the same spirit 
have been developed that incorporate heating due to accretion 111221 . but the full 
problem requires numerical treatment. Several codes are available for its efficient 
solution 1551 ITOOl 12941 13151 . 

If our main interest is in the physical conditions at the disk mid-plane within the 
normal planet forming region (i.e. excluding very large orbital distances where the 
disk is becoming optically thin), or in the continuum Spectral Energy Distribution 
(SED), models for the dust temperature largely suffice. Dust dominates both the 
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absorption of starlight and thermal emission of reprocessed radiation and accretion 
heating, and in the mid-plane the gas and dust temperatures are normally closely 
equal (to within ~ 10% for an optical extinction Ay >0.1 Il68l l. The gas tempera¬ 
ture near the surface of the disk is, however, of critical importance for a number of 
applications, 

• Interpretation of sub-mm data, where the observable emission is rotational transi¬ 
tions of molecules such as CO and HCO + . These observations frequently probe 
the outer regions of disks, at depths where the molecules are not (of course!) 
photo-dissociated but where the gas is warm and not in equilibrium with the 
dust. 

• Interpretation of near-IR and far-IR data, often from the inner disk, where we are 
seeing ro-vibrational transitions of molecules along with fine-structure cooling 
lines such as [CII] and [OI]. 

• Chemistry. It’s cold at the disk mid-plane, and chemical reactions are sluggish. 
Although the densities are much lower in the disk atmosphere, the increased tem¬ 
peratures and exposure to higher energy stellar photons make the upper regions 
of the disk important for chemistry 111 441 . 

The properties of gas near the surface of disks are very closely tied to the incident 
flux of ultraviolet radiation from the star. Stellar UV radiation ionizes and dissoci¬ 
ates atoms and molecules, and heats the gas by ejecting electrons from dust grains 
(grain photoelectric heating). Depending upon the temperature and density, the heat¬ 
ing is balanced by cooling from rotational transitions of molecules (especially CO) 
and atomic fine structure lines. Also important is energy exchange due to inelastic 
collisions between gas molecules and dust particles (thermal accommodation) — if 
this process is too efficient the gas temperature will revert to match the dust which 
is absorbing and emitting the bulk of the star’s bolometric luminosity. 

Photoelectric heating ll94l is typically the dominant process for dusty gas ex¬ 
posed to an ultraviolet radiation field. The work function of graphite grains (the 
minimum energy required to free an electron from them) is around 5 eV, so 10 eV 
FUV photons can eject electrons from uncharged grains with 5 eV of kinetic en¬ 
ergy that ultimately heats the gas. Ejection occurs with a probability of the order of 
0.1, so the overall efficiency (the fraction of the incident FUV energy that goes into 
heating the gas rather than the dust) can be rather high, around 5%. 

A detailed evaluation of the photoelectric heating rate is involved, and resistant 
to a fully first-principles calculation. Weingartner & Draine 13481 give a detailed 
description. Here, we sketch the main principles following Kamp & van Zadelhoff 
11691 . who developed models for the gas temperature in A star disks. We consider 
a stellar radiation spectrum F v impinging on grains of graphite (work function w = 
4.4 eV 113481 ) and silicate (w = 8.0 ev). For micron-sized grains the work function, 
which is a property of the bulk material, is equivalent to the ionization potential — 
the energy difference between infinity and the highest occupied energy level in the 
solid. Additionally, the probability for absorption when a photon strikes a grain is 
2abs ~ 1 • Most absorbed photons, however, do not eject electrons, rather their energy 
goes entirely into heating the dust grain. The yield of emitted electrons is some 
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function of photon energy Y ( hv ), and they have some spectrum of kinetic energy E, 
roughly described by |94l f(E,hv) « (hv — w) _1 . If the grains are charged (e.g. by 
prior emission of photoelectrons) then the kinetic energy (E — eU) available to heat 
the gas is that left over once the electron has escaped the electrostatic potential eU 
of the grain. The heating rate is then Il69l . 



(31) 


where n/j is the number density of hydrogen atoms, cr is the geometric cross-section 
per hydrogen nucleus, and the lower limits express the minimum frequency v,h of 
a photon that can overcome the work function and the minimum energy of a pho¬ 
toelectron that can escape from a charged grain. An assessment of the photoelec¬ 
tric heating rate then requires knowledge of the functions Y and /, specification 
of the radiation field F v , and calculation of the typical charge on grains of differ¬ 
ent sizes 1%113481 . The physics is conceptually identical but quantitatively distinct 
when the grains in question are extremely small (e.g. Polycyclic Aromatic Hydro¬ 
carbons, PAHs) Il35ll348l . 

The rate of energy exchange from inelastic gas-grain collisions can be calculated 
with a collision rate argument. Consider grains with geometric cross-section Oj = 
7r(s 2 ) and number density n ( \, colliding with hydrogen atoms with number density 
«//. The thermal speed of the hydrogen atoms is v* = (8 ksTg/nmnY^ 2 and the 
average kinetic energy of the molecules on striking the surface is 2ksT g . The cooling 
rate per unit volume due to gas-grain collisions can then be written in the form ll64| . 



(32) 


where T g and 7)/ are the temperatures of the gas and dust respectively. The subtleties 
of the calculation are reflected in the “accommodation co-efficient” ttr , which is 
typically O.j ss 0.3 for silicate and carbon grains. For a specified volumetric heating 
rate (and assumptions as to the gas to dust ratio and properties of the grains), this 
expression can be used to estimate the density below which the thermal properties 
of gas and dust decouple. 

In addition to cooling that occurs indirectly, as a consequence of gas-grain colli¬ 
sions, gas in the upper layers of disks also cools radiatively. In the molecular layer of 
the disk, the dominant coolant is typically CO, as this is the most abundant molecule 
that is not homonuclear (diatomic molecules, such as Hi, have no permanent elec¬ 
tric dipole moment and hence radiate inefficiently). At higher temperatures — only 
attained in the very rarefied uppermost regions of the disk atmosphere — cooling by 
Lya emission becomes important. Qualitatively, there are then three distinct layers 
in the disk: 

• A cool mid-plane region, where dust and gas have the same temperature and dust 

cooling is dominant. 
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Fig. 5 Illustration of some of the physical processes determining the temperature and emission 
properties of irradiated protoplanetary disks. 


• A warm surface layer, in which both dust and gas have temperatures that exceed 
the mid-plane value. The gas in the warm layer can be substantially hotter than 
the dust (T ~ 10 3 K at 1 AU), and cools both by dust-gas collisions and by CO 
rotational-vibrational transitions. 

• A hot, low-density atmosphere, where Lya radiation and other atomic lines (e.g. 
0[I]) cool the gas. 

The disk structure that results from these heating and cooling processes is illustrated 
in Figure [5] 


2.3 Ionization structure 

The degree of ionization of the gas in protoplanetary disks is important because 
it is key to understanding how gas couples to magnetic fields, and thence to un¬ 
derstanding the role of magnetic fields in the formation of disks, in the sustenance 
of turbulence within them, and in the generation of jets and magnetohydrodynamic 
(MHD) winds. At the most basic level, we care about the ratio of the number density 
of free electrons n e to the number density of neutrals. 
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n„ 

though we should remember that dust grains can also bear charges and carry cur¬ 
rents. We will consider separately the thermodynamic equilibrium process of ther¬ 
mal (or collisional) ionization, which typically dominates above T ~ 10 3 K, and 
non-thermal ionization due to photons or particles that have an energy well in ex¬ 
cess of the typical thermal energy in the gas. 

In anticipation of results that will be derived in £j4} we note that very low and 
seemingly negligible levels of ionization — x e 10 10 - often suffice to couple 
magnetic fields to the fluid. We need to worry about small effects when considering 
ionization. 


2.3.1 Thermal ionization 


Thermal ionization of the alkali metals is important in the innermost regions of the 
disk, usually well inside 1 AU. In thermal equilibrium the ionization state of a single 
species with ionization potential % is obeys the Saha equation 12991 . 


2 U ion ( 2nm e k B T 

u v ^ 


3/2 

exp[-*Afl7l. 


(34) 


Here, « lon and n are the number densities of the ionized and neutral species, and 
n e (= n lon ) is the electron number density. The partition functions for the ions and 
neutrals are U' on and U, and the electron mass is m e . The temperature dependence 
is not quite just the normal exponential Boltzmann factor, because the ionized state 
is favored on entropy grounds over the neutral state. 

In protoplanetary disks thermal ionization becomes significant when the tem¬ 
perature becomes high enough to start ionizing alkali metals. For potassium, the 
ionization potential % = 4.34 eV. We write the abundance of potassium relative to 
all other neutral species as / = riK/n n , and define the ionization fraction x, 

x =—. (35) 

n„ 


While potassium remains weakly ionized, the Saha equation gives, 


x — 1(T 


12 



- 1/2 / t \ 3/4 exp[—2.52 x 10 4 /T] 
V10 3 kJ 1.14 x 10- 11 


(36) 


where the final numerical factor in the denominator is the value of the exponent at 
10 3 K. The ionization fraction at different temperatures is shown in Figure^ Ioniza¬ 
tion fractions large enough to be interesting for studies of magnetic field coupling 
are reached at temperatures of T ~ 10 3 K although the numbers remain extremely 
small - of the order of x ~ 10 13 for these parameters. 
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Fig. 6 The thermal ionization fraction as a function of temperature predicted by the Saha equation 
for the inner disk. Here we assume that potassium, with ionization potential % = 4.34 eV and 
fractional abundance / = 10 7 , is the only element of interest for the ionization. The number 
density of neutrals is taken to be n„ = 10 15 cm~ 3 . 


2.3.2 Non-thermal ionization 

Outside the region close to the star where thermal ionization is possible, any rem¬ 
nant levels of ionization are controlled by non-thermal processes. Considerations of 
thermodynamic equilibrium are not relevant, and we need to explicitly balance the 
rate of ionization by high-energy particles or photons against the rate of recombina¬ 
tion within the disk gas. 

There is no shortage of potentially important sources of ionization. Ordering 
them roughly in order of their penetrating power, ideas that have been suggested 
include, 

• Ultraviolet photons (from the star, or from other stars in a cluster) 

• Stellar X-rays 

• Cosmic rays 

• Energetic protons from a stellar corona 113338 

• Particles produced from radioactive decay of nuclides within the disk B316il 

• Electric discharges. B240B 

We will limit our discussion to the first three of these processes. 

The coronae of T Tauri stars are powerful sources of keV X-rays t273l . Typical 
luminosities are Lx — 10 28 — 10 31 erg s 1 , in X-rays with temperatures ksTx of a 
few keV. The physics of the interaction of these X-rays with the disk gas involves 
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Compton scattering and absorption by photo-ionization, which has a cross-section 
<7 ~ 10 22 cm 2 for keV energies, decreasing with photon energy roughly as L ph l r 
Given an input stellar spectrum and assumptions as to where the X-rays originate, 
the scattering and absorption physics can be calculated using radiative transfer codes 
to deduce the ionization rate within the disk llOll . 

Depending upon the level of detail needed for a particular application, the re¬ 
sults of numerical radiative transfer calculations can be approximated analytically. 
For a relatively hard stellar spectrum (kgTx = 5 keV), the ionization rate fairly deep 
within the disk scales with radius r and vertical column from the disk surface E 
as r“ 2 exp[— E/S g cm' 2 ] 13341 . A more detailed fit is given by Bai & Goodman 
(2009) li29l . For an X-ray luminosity scaled to Lx ,29 = Lx/10 29 erg s' 1 they repre¬ 
sent the numerical results with two components, 

L^ (liu)^ 2 = Clex PHW)°] + beM-(E/E 2 ) p ] + ... (37) 

where E is the vertical column density from the top of the disk and symmetric 
terms in the column density from the bottom of the disk are implied. For ky{I\ = 
3 keV and Solar composition gas the fit parameters are ^ = 6 x 10' 12 s' 1 , Ci = 
10' 15 s' 1 , El = 3.4 x 10' 3 g cm' 2 , E 2 = 1.59 g cm' 2 , a = 0.4 and J3 = 0.65. 
For kgTx = 5 keV the fit parameters are = 4 x 10' 12 s' 1 , ^=2x 10' 15 s' 1 , 
E\ = 6.8 x 10' 3 g cm' 2 , E 2 = 2.27 g cm' 2 , a = 0.5 and p = 0.7. 

Figure FT] shows the estimates for £(£) for a stellar X-ray luminosity of Lx = 
10 30 erg s incident on the disk at 1 AU. If one is mainly interested in regions 
of the disk more than « 10 g cm' 2 away from the surfaces the single exponential 
fit given by Turner & Sano (2008) j334l may suffice. The more complex fitting 
function given by equation ( [37] ) captures the much higher rates of ionization due to 
X-rays higher up in the disk atmosphere. 

Cosmic rays are another potential source of disk ionization. A standard descrip¬ 
tion of the interstellar cosmic ray flux gives them an unattenuated ionization rate of 
C cr ~ 10' 17 — 10 16 s' 1 and an exponential stopping length of 96 g cm' 2 (substan¬ 
tially greater than even high energy stellar X-raysWith these parameters. X-rays 
would remain the primary source of ionization in the upper ss 50 g cm' 2 of the disk, 
but cosmic rays would dominate in the region between about 50 and 500 g cm' 2 . 
It is unclear, however, whether the unattenuated interstellar medium flux of cos¬ 
mic rays typically reaches the surfaces of protoplanatary disks. The magnetic fields 
embedded in the Solar wind form a partial barrier to incoming cosmic ray parti¬ 
cles, whose effect is seen in a modulation of the observed flux with the Solar cycle. 
T Tauri stars could have much stronger stellar winds that exclude cosmic rays effi¬ 
ciently. Indeed, chemical modeling of molecular line data suggests that cosmic rays 
are substantially excluded (to a level Ccr ~ 10 I9 ) from the disk around the nearby 
star TW Hya lf83l . though how pervasive this phenomenon is remains unknown. If 


- Umebayashi & Nakano noted that if cosmic rays have an approximately isotropic angular distri¬ 
bution at the disk surface, geometric effects lead to a faster than exponential attenuation deep in 
the disk f336l . 
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Fig. 7 Estimates of the non-thermal ionization rate due to X-rays (red curves), unshielded cosmic 
rays (blue) and radioactive decay of short-lived nuclides (blue), plotted as a function of the vertical 
column density from the disk surface. The solid red curve shows the Bai & Goodman (2009) result 
for an X-ray temperature kgTx = 5 keV, the dashed curve their result for kgTx = 3 keV. The dot- 
dashed red curve shows a simpler formula proposed by Turner & Sano (2008). All of the X-ray 
results have been normalized to a flux of Lx = 10 30 erg s _1 and a radius of 1 AU. 


cosmic rays are not present, the only guaranteed source of ionization at columns 
more than « 100 g cm~ 2 away from the disk surfaces is radioactive decay. 

If our main interest is in conditions at r ~ 1 AU the surface density in gas is 
typically E ~ 10 3 gcm~ 2 and X-rays, which are our main concern, will not reach the 
mid-plane. The situation is different further out. At 100 AU typical surface densities 
are much lower — 1 g cm -2 might be reasonable — and X-rays will sustain a non¬ 
zero rate of ionizations throughout that column. On these scales ultraviolet photons 
can also be important. Stellar FUV radiation will ionize carbon and sulphur atoms 
near the disk surface, yielding a relatively high electron fraction x e ~ 10 5 . The 
ionized skin that results is shallow, penetrating to a vertical column of just 0.01 — 
0.1 g cm -2 , but enough to be significant in the tenuous outer disk 12661 . 

Whether cosmic rays are excluded from most or all T Tauri disks is an im¬ 
portant question which may be addressed observationally in the near future. 
Theoretically one should also remember that the amount of power involved 
in non-thermal ionization is rather small when compared to that liberated by 
accretion 1)1561 . Any internal disk process that could convert a small fraction 









Physical Processes in Protoplanetary Disks 


19 


of the accretion energy into non-thermal particles would likely matter for the 
ionization state. 


As with ionization, the rate of recombination within the disk can be calculated 
from complex numerical models that track reactions (often numbering in the thou¬ 
sands) between dozens of different species. The following discussion, which bor¬ 
rows heavily from the description given by Ilgner & Nelson (2006) EH and Fro- 
mang (2013) (TT3l . is intended only to outline some of the important principles. At 
the broadest level of discussion we need to consider gas-phase recombination reac¬ 
tions (involving molecular and gas-phase metal ions) along with recombination on 
the surface of dust grains. 

The principles of gas-phase recombination can be illustrated by considering the 
possible reactions between electrons and generic molecules m and metal atoms M 
1 25211152 1. The basic reactions are then, 

• Ionization, 

m—>-m + + e _ , (38) 

with rate £. A specific example is kb —> Hj +e~. 

• Recombination with molecular ions, 

m + +e~—>■ m, (39) 

with rate a = 3 x 10 _6 !T _1 / 2 cm 3 s ' 1 . An example is the dissociative recombi¬ 
nation reaction HCO + + e —> CO + H. 

• Recombination with gas-phase metal ions, 

M + +e” — M + hv, (40) 

with rate y = 3 x 10 ~ 11 7 ’“ 1 / 2 cm 3 s -1 . An example is Mg + + e~ —> Mg + hv. 

• Charge exchange reactions. 


m + + M-S>m + M + , (41) 

with rate j3 = 3 x 10 9 cm 3 s 1 . An example is HCO + + Mg —» Mg + + HCO. 

From such a set of reactions we form differential equations describing the time 
evolution of the number density of species involved. For the molecular abundance 
/ 7 m , for example, we have, 

d fi 

-7— = -C /7 m + Ctn e n m + +priMn m + = 0, (42) 

at 

where the second equality follows from assuming that the system has reached equi¬ 
librium. The resulting system of algebraic equations has simple limiting solutions. 
For example, if there are no significant reactions involving metals then the above 
equation, together with the condition of charge neutrality (n m + = n e ), gives an elec- 
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tron fraction x e = n e /n m . 



(43) 


In the more general case, the network yields a cubic equation which can be solved 
for the electron fraction as a function of the gas-phase metal abundance HSU. Typ¬ 
ically the presence of metal atoms and ions is important for the ionization level. 

Recombination can also occur on the surfaces of dusty or icy grains. The simplest 
reactions we might consider are. 


e + gr —► gr 

m + + gr~ —► gr + m. (44) 

If the first of these reactions is rate-limiting, then we can write a modified version 
of the ordinary differential equation (equation |42] > that includes grain processes. 
Ignoring metals for simplicity, 

dr/ in , . . c . 

=-Qn m + an e n m + + Ov e n gT n e ~, (45) 

where cr is the cross-section of grains to adhesive collisions with free electrons 
and v e is the electron thermal velocity. In the limit where only grains contribute to 
recombination we then find, 

Xe = ---■ (46) 

If the grains are mono-disperse with radius s, then x e 1 /On g: oc s, and recombina¬ 
tion on grain surfaces will be more important for small grain sizes. We also note that 
the dependence on the ionization rate £ is linear, rather than the square root depen¬ 
dence found in the gas-phase case. As for metals, grain populations with commonly 
assumed size distributions are found to matter for the ionization level. 

The above discussion of recombination leaves a great deal unsaid. For grains, an 
important additional consideration is related to the typical charge state, which needs 
to be calculated (96). A good comparison of different networks for the calculation of 
the ionization state is given by Ilgner & Nelson (2006) II152H . while Bai & Goodman 
(2009) li29l provide a clear discussion of the important processes. 

The accuracy needed from calculations of ionization equilibrium is strongly 
problem-dependent (and in some cases, such as if we don’t know if cosmic 
rays are present for a particular system, high accuracy may be illusory). The 
reader who encounters a problem involving the ionization level is advised 
to consider whether a simple analytic approximation is adequate, or whether 
solution of a full chemical network is required. Both new analytic models 
(that include dust), and flexible freely-available codes for computing detailed 
solutions, would be valuable. 
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3 Disk evolution 

The population of protoplanetary disks is observed to evolve, but why this should 
happen remains hard to fully explain. For a geometrically thin, low-mass disk, the 
deviation from a point-mass Keplerian rotation curve is small (c.f. equation 1 13 [> and 
the specific angular momentum. 


l{r) = r 2 £ 2k = \JGM*r « r 1 / 2 


(47) 


is an increasing function of orbital radius. To accrete, gas in the disk must lose angu¬ 
lar momentum, and the central theoretical problem in disk evolution is to understand 
this process. 

Within any shearing fluid, momentum is transported in the cross-stream direction 
because the random motion of molecules leads to collisions between particles that 
have different velocities. The classical approach to disk evolution 8219112761 treats 
the disk as a vertically thin axisymmetric sheet of viscous fluid, and leads to a fairly 
simple equation for the time evolution of the disk surface density E(r,t). There 
appears to be a fatal flaw to this approach, because the molecular viscosity of the 
gas is much too small to lead to any significant rate of disk evolution. But it’s not 
as bad as it seems. The classical disk evolution equation involves few assumptions 
beyond the immutable laws of mass and angular momentum conservation, and as we 
shall see is therefore approximately valid if the “viscosity” is re-interpreted as the 
outcome of a turbulent process. We will have (much) more to say about the possible 
origin of disk turbulence in (J4] 

Redistribution of angular momentum within the gas disk is not the only route 
to evolution. An almost equally long-studied suggestion [57l| is that gas accretes 
because a magnetohydrodynamic (MHD) wind removes angular momentum entirely 
from the disk. Winds and viscosity have frequently between seen as orthogonal and 
competing hypotheses for disk evolution, but there is evidence suggesting that both 
processes are simultaneously important in regions of protoplanetary disks. 


3.1 The classical equations 

The evolution of a flat, circular and geometrically thin (( h/r ) <C 1) viscous disk 
follows from the equations of mass and angular momentum conservation 82768 . 
Given a surface density E(r, t ), radial velocity v r (r, t) and angular velocity Q (r), the 
continuity equation in cylindrical co-ordinates yields. 



( 48 ) 


Angular momentum conservation gives. 
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r^- (r 2 QE) + (r 2 £2 ■ rEv r 

dt dr 


1 dG 
2n dr ' 


(49) 


where £2 (r) is time-independent but need not be the point mass Keplerian angular 
velocity. The rate of change of disk angular momentum is given by the change in 
surface density due to radial flows and by the difference in the torque exerted on the 
annulus by stresses at the inner and outer edges. For a viscous fluid the torque G has 
the form, 

G = 2nr-vEr-—r (50) 

dr 


where v is the kinematic viscosity. The torque is the product of the circumference, 
the viscous force per unit length, and the lever arm r, and scales with the gradient 
of the angular velocity. 

To obtain the surface density evolution equation in its usual form we first elimi¬ 
nate v r by substituting for dE/dt in equation (49 1 from equation (48 i. This gives an 
expression for rEv r , which we substitute back into equation (|48]> to yield. 


dE 

dt 


1 d 

r dr 


1 d 

( r 2 £2)' dr 


(vEr'Q 1 ) 


(51) 


where the primes denote differentiation with respect to radius. Specialize to a point 
mass Keplerian potential (£2 r 3/2 ) we then find that viscous redistribution of 
angular momentum within a thin disk obeys an equation. 


dE 

dt 


3 d 
r dr 



(52) 


This equation is a diffusive partial differential equation for the evolution of the gas, 
which has a radial velocity. 


IV 


3 d 
Er^-I- dr 



(53) 


The equation is linear if the viscosity v is independent of E. 

Some useful rules of thumb for the rate of evolution implied by equation (52 1 can 
be deduced with a change of variables. Defining, 


X = 2 r 1 / 2 

(54) 

/ = -£x, 

j 2 

(55) 


and taking the viscosity v to be constant, we get a simpler looking diffusion equa¬ 
tion, 


df 

dt 


= D 


*1 

dx 2 


(56) 


with a diffusion coefficient D, 


D = 


12v 

W 


(57) 
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Fig. 8 The time-dependent solution to the disk evolution equation with v = constant, showing the 
spreading of a ring of gas initially orbiting at r = tq. From top down the curves show the surface 
density as a function of the scaled time variable T = Hvr^ 2 !, for T = 0.004, T = 0.008, T = 0.016, 
% = 0.032, t = 0.064, t = 0.128, and z = 0.256. 


The diffusion time scale across a scale AX for an equation of the form of equa¬ 
tion (561 is just ( AX) 1 /D. Converting back to the physical variables, the time 
scale on which viscosity will smooth surface density gradients on a scale Ar is 
T v ~ (Ar) 1 /v. For a disk with characteristic size r, the surface density at all radii 
will evolve on a time scale. 



(58) 


This is the viscous time scale of the disk. 

We can gain some intuition into how equation ( [52] ) works by inspecting time- 
dependent analytic solutions that can be derived for special forms of the viscosity 
v(r). For v = constant a Green’s function solution is possible. Suppose that at t = 0 
the gas lies in a thin ring of mass m at radius rp. 


£(r,t = 0) = 


m 

2nro 


8(r 


ro), 


(59) 


where 5(r —rp) is a Dirac delta function. With boundary conditions that impose 
zero-torque at r = 0 and allow for free expansion toward r = °° the solution is 12191 . 


£(x,t) 


m 1 
Kr^ T 


exp 


(1+x 2 )' 


^ 1/4 



(60) 












24 


Philip J. Armitage 


in terms of dimensionless variables x = r/ro, T = 12vr () 2 t, and where I\u is a mod¬ 
ified Bessel function of the first kind. This solution is plotted in Figure [8] and il¬ 
lustrates generic features of viscous disk evolution. As t increases the ring spreads 
diffusively, with the mass flowing toward r = 0 while the angular momentum is car¬ 
ried by a negligible fraction of the mass toward r = °°. This segregation of mass and 
angular momentum is generic to the evolution of a viscous disk, and must occur if 
accretion is to proceed without overall angular momentum loss. 


3.1.1 Limits of validity 


Protoplanetary disks are not viscous fluids in the same way that honey is a viscous 
fluid (or, for that matter, in the same way as the mantle of the Earth is viscous). 
To order of magnitude precision, the viscosity of a gas v ~ v^A, where i', h is the 
thermal speed of the molecules and the mean-free path A is. 


1 

no 


(61) 


Here n is the number density of molecules with collision cross-section o. Taking 
cr to be roughly the physical size of a hydrogen molecule, o ~ k( 10 8 cm) 2 , and 
conditions appropriate to 1 AU ( n ~ 10 15 cm~ 3 , v t h ~ 10 5 cm s -1 ) we estimate. 


A ~ 3 cm 

v ~ 3 x 10 5 cm 2 s *. 


(62) 


This is not a large viscosity. The implied viscous time according to equation (58 i 
is of the order of 10 13 yr, far in excess of observationally inferred time scales of 
protoplanetary disk evolution. If we nevertheless press on and use equation ( [52| to 
model disk evolution, we are implicitly modeling a system that is not an ordinary 
viscous fluid with a viscous equation. We need to understand when this is a valid 
approximation. 

The first possibility, introduced by Shakura & Sunyaev (1973) in their paper on 
black hole accretion disks 13031 . retains the idea that angular momentum is con¬ 
served within the disk system, but supposes that turbulence rather than molecular 
processes is the agent of angular momentum transport. Looking back at the deriva¬ 


tion of the disk evolution equation (52 1 , we note that the fluid properties of molec¬ 
ular viscosity only enter twice, (i) in the specific expression for G (which, e.g. is 
linear in the rate of shear) and (ii) in the more basic assumption that angular mo¬ 
mentum transport is determined by the local fluid properties. The rest of the deriva¬ 
tion involves only conservation laws that hold irrespective of the nature of transport. 
Plausibly then, a disk in which angular momentum is redistributed by the action 
of turbulence should still be describable by a diffusive equation, provided that the 
turbulence is a local process. Proceeding rigorously, Balbus & Papaloizou (1999) 
140) showed that MHD turbulence is in principle local in this sense, whereas angu¬ 
lar momentum transport by self-gravity is in principle not. At the level of the basic 
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axisymmetric evolution equation then — before we turn to questions of what deter¬ 
mines v, or how boundary layers behave where the shear is reversed — we have not 
committed any cardinal sin in starting from the viscous disk equation. 

Greater care is needed in situations where the disk flow is no longer axisym¬ 
metric. Fluids obey the Navier-Stokes equations, but there is no guarantee that a 
turbulent disk with a complex geometry will behave in the same way as a viscous 
Navier-Stokes flow with effective kinematic and bulk viscosities. In eccentric disks, 
for example, even the most basic properties (such as whether the eccentricity grows 
or decays) depend upon the nature of the angular momentum transport 112481 . 

The disk evolution equation will also need modification if there are external 
sources or sinks of mass or angular momentum. If the disk gains or loses mass 
at a rate E ( r , t), and if that gas has the same specific angular momentum as the disk , 
then the modification is trivial. 



(63) 


Disk evolution in the presence of thermally driven winds (such as photo-evaporative 
flows) can be described with this equation. Alternatively, we may consider a disk 
subject to an external torque that drives a radial flow with velocity v r . ex i. This adds 
an advective term. 



(64) 


The qualitative evolution of the disk — for example the tendency for the outer re¬ 
gions to expand to conserve angular momentum — can be changed if there is even 
a modest external torque on the system. 

The evolution of disk populations is frequently modeled using the diffusion 
equation p2| ). In addition to the formal question of whether this is mathe¬ 
matically sensible, it is not entirely obvious that protoplanetary disks evolve 
significantly under viscous torques at all radii. In the outer regions, especially, 
it is possible that the initial surface density distribution is modified more by 
thermal or magnetic winds. 


3.1.2 The a prescription 

Molecular viscosity depends in a calculable way upon the density, temperature and 
composition of the fluid. Can anything similar be said about the “effective viscosity” 
present in disks? The standard approach is to write the viscosity as the product of 
the characteristic velocity and spatial scales in the disk. 


v = acji, 


(65) 
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where a is a dimensionless parameter. This ansatz (introduced in a related form in 
13031 ) is known as the Shakura-Sunyaev a prescription. 

We can view the a prescription in two ways. The “weak” version is to regard it 
as a re-parameterization of the viscosity that (perhaps) describes the leading order 
scaling expected in disks (so that a is a more slowly varying function of tempera¬ 
ture, radius etc than v). This is useful, and along with convention is the reason why 
numerical simulations of turbulent transport are invariably reported in terms of an 
effective a. One can also adopt a “strong” version of the prescription in which a 
is assumed to be a constant. This is powerful, as it allows for the development of 
a predictive theory of disk structure that is based on only one free parameter (for a 
textbook discussion see Frank, King & Raine cm, or for an application to proto¬ 
planetary disks see, e.g. 109]). However, its use must be justified on a case by case 
basis, as there is no reason why a should be a constant. Constant a models proba¬ 
bly work better in highly ionized disk systems, where angular momentum transport 
across a broad range of radii occurs via the simplest version of magnetorotational 
instability (38], than in protoplanetary disks where the physical origin of angular 
momentum transport is more complex m. 


3.2 Boundary conditions 

Solving equation ( [52] ) requires the imposition of boundary conditions. The most 
common, and simplest, is a zero-torque inner boundary condition, which exactly 
conserves the initial angular momentum content of the disk. If the star has a dynam¬ 
ically significant magnetic field, however, or if the disk is part of a binary system, 
other boundary conditions may be more appropriate. 


3.2.1 Zero-torque boundary conditions 


A steady-state solution to equation ( |52| ) with a zero-torque inner boundary condition 
is derived by starting from the angular momentum conservation equation (Equa¬ 
tion F 


491. Setting the time derivative to zero and integrating we have. 


2nrT.v r ■ r z £2 


: 2nr vL - \- constant. 

dr 


( 66 ) 


In terms of the mass accretion rate M = —2nrEv r we can write this in the form, 

. n "3 dd 

— M ■ r z Q = 2nrvT. -bconstant, (67) 

dr 

where the constant of integration, which is an angular momentum flux, is as yet 
undermined. To specify the constant, we note that if there is a point where df2 /dr = 
0 the viscous stress vanishes, and the constant is just the advective flux of angular 
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r bl -dQ / dr = 0 



r = R, 


Fig. 9 A sketch of what the angular velocity profile £2 ( r ) must look like if the disk extends down 
to the surface of a slowly rotating star. By continuity there must be a point — usually close to the 
stellar surface — where /dr = 0 and the viscous stress vanishes. 


momentum advected, 

constant = —M ■ r 2 £2. (68) 

The physical situation where A£2/Ar = 0 is where the protoplanetary disk extends 
all the way down to the surface of a slowly rotating star. The disk and the star 
form a single fluid system, and the angular velocity (shown in Figure |9]> must be a 
continuous function that connects £2 = 0 in the star to £2 r 3 / 2 within the disk. The 
viscous stress must then vanish at some radius R, + /'hi, where ru is the width of the 
boundary layer that separates the star from the Keplerian part of the disk. Within the 
boundary layer the angular velocity increases with radius, and gravity is balanced 
against a combination of rotation and radial pressure support. Elementary arguments 
12751 show that in many case the boundary layer is narrow, so that R., + ~ R,. 

We then find that, 

■ , /gm* 

constant ~ -MR 2 J —(69) 
and the steady-state solution for the disk simplifies to, 



For a specified viscosity this equation gives the steady state surface density profile of 
a disk with a constant accretion rate M. Away from the inner boundary E(rj ^ v 1 , 
and the radial velocity (Equation [53} is v r = —3v/2r. 
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In obtaining equation (70 1 we have derived an expression for a Keplerian disk 
via an argument that relies on the non-Keplerian form of Q (r) in a boundary layer. 
The resulting expression for the surface density is valid in the disk at r^> R.,, but 
would not work well close to the star even if there is a boundary layer. To model the 
boundary layer properly, we would need equations that self-consistently determine 
the angular velocity along with the surface density 1272ft . 


3.2.2 Magnetospheric accretion 

For protoplanetary disks the stellar magnetic field can have a dominant influence on 
the disk close to the star 111781 . The simplest magnetic geometry involves a dipolar 
stellar magnetic field that is aligned with the stellar rotation axis and perpendicular 
to the disk plane. The unperturbed field then has a vertical component at the disk 
surface, 



In the presence of a disk,the vertical field will thread the disk gas be distorted by 
differential rotation between the Keplerian disk and the star. The differential rotation 
twists the field lines that couple the disk to the star, generating an azimuthal field 
component at the disk surface B () and a magnetic torque per unit area (counting both 
upper and lower disk surfaces). 


T = 


B z B<j) 

~2n 


(72) 


Computing the perturbed field accurately is hard (for simulation results see, e.g. 
12951 ). but it is easy to identify the qualitative effect that it has on the disk. Defining 
the co-rotation radius r co as the radius where the field lines have the same angular 
velocity as that of Keplerian gas in the disk. 


r CO 


( gm*p 2 \ 1/3 

V 47 1 2 J 


(73) 


there are two regions of star-disk magnetic interaction: 

• Interior to co-rotation (r < r co ) the disk gas has a greater angular velocity than 
the field lines. Field lines that link the disk and the star here are dragged forward 
by the disk, and exert a braking torque that removes angular momentum from the 
disk gas. 

• Outside co-rotation (r > r co ) the disk gas has a smaller angular velocity than 
the field lines. The field lines are dragged backward by the disk, and there is a 
positive torque on the disk gas. 

Young stars are typically rapid rotators |f62l , so the co-rotation radius lies in the 
inner disk. For P = 7 days, for example, the co-rotation radius around a Solar mass 
star is at r co ~ 15 Rq or 0.07 AU. 
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The presence of a stellar magnetic torque violates the assumption of a zero-torque 
boundary condition, though the steady-state solution we derived previously (equa¬ 
tion [70]) will generally still apply at sufficiently large radius. The strong radial de¬ 
pendence of the stellar magnetic torque means that there is only a narrow window of 
parameters where the torque will be significant yet still allow the disk to extend to 
the stellar surface. More commonly, a dynamically significant stellar field will dis¬ 
rupt the inner disk entirely, yielding a magnetospheric regime of accretion in which 
the terminal phase of accretion is along stellar magnetic held lines. The disruption 
(or magnetospheric) radius r m can be estimated in various ways 111781 . but all yield 
the same scaling as the spherical Alfven radius. 


\ 2/? 

My/GM~J ’ 


(74) 


where B* is the stellar surface held (dehned such that B, is the dipole moment) 
and k a constant of the order of unity. Taking k = 1 for a Solar mass star. 
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(75) 


Often, the magnetospheric radius is comparable to the co-rotation radius. 


3.2.3 Accretion on to and in binaries 


Boundary conditions for disk evolution also need modihcation in binary systems. 
For a coplanar disk orbiting interior to a prograde stellar binary companion, tidal 
torques from the companion remove angular momentum from the outer disk and 
prevent it from expanding too far 126211 . The tidal truncation radius roughly corre¬ 
sponds to the largest simple periodic orbit in the binary potential 112601 . which is at 
about 40% of the orbital separation for a binary with mass ratio Mi/M\ = 0.5. The 
tidal torque is a function of radius, but to a hrst approximation one may assume that 
tides impose a rigid no-expansion condition at r = r out . From equation (|53j), 


d 

dr 



r=r 0 ut 


= 0 . 


(76) 


This type of boundary condition may also be an appropriate approximation for a 
circumplanetary disk truncated by stellar tides 12291 . 

The size of the disk (and even whether it is tidally truncated at all) will be 
different if the disk is substantially misaligned with respect to the orbital plane 
of the binary 1214. 237 1. 
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Fig. 10 The time-dependent analytic solution (2771 to the disk evolution equation with a v r = 0 
boundary condition at r = 1 for the case v = r. The curves show the evolution of a ring of gas 
initially at r = 2, at times f = 0.002, t = 0.004, t = 0.008 etc. The bold curve is at t = 0.128, and 
dashed curves show later times. Gas initially accretes, but eventually decretes due to the torque 
being applied at the boundary. 


An exterior circumbinary disk will also experience stellar gravitational torques, 
which in this case add angular momentum to the disk and slow viscous inflow. How 
best to model these torques is an open question, particularly in the case of extreme 
mass ratio binaries composed of a star and a massive planet. Pringle (1991) derived 
an illuminating analytic solution for circumbinary disk evolution izm under the 
assumption that tidal torques completely prevent inflow past some radius r = r m . 
With this assumption the boundary condition at r = r ln is v r = 0, and the task is to 
find a solution to equation ( [52] ) with this finite radius boundary condition. A simple 
solution is possible for v = kr, with k a constant. Defining scaled variables, 

x = r 1 / 2 

(7 = Zr 3 / 2 , (77) 

the t > 0 solution for an initial mass distribution, 

cr(x,f = 0) = Oo<5(a — -71), (78) 


is 12771 . 
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(79) 


a = 
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The solution, plotted in Figure[lO] can be compared to the zero-torque solution (Fig- 
ure[8] though note this is for a constant viscosity). The initial evolution is similar, but 
at late times the torque which precludes inflow past r m causes qualitatively different 
behavior. The disk switches from an accretion to a decretion disk, with an outward 
flow of mass driven by the binary torque. 

Disks resembling the classical decretion disk solution may be realized under 
some circumstances, for example around rapidly rotating and strongly magnetized 
stars. Around binaries, however, it appears hard to completely shut off inflow. Nu¬ 
merical simulations show that angular momentum transfer to the disk co-exists with 
persistent inflow into a low density cavity containing the binary l20j[93 l. How best 
to represent this complexity in a one-dimensional model is not entirely obvious. 


3.3 Viscous heating 


Although stellar irradiation is often the dominant source of heat for protoplanetary 
disks ( §2.2| l, dissipation of gravitational potential energy associated with accretion 
is also important. Ignoring irradiation for the time being, we can derive the effective 
temperature profile of a steady-state viscous disk. In the regime where the classical 
equations are valid, the fluid dissipation per unit area is 12761 . 


Q + = 9 -v£L l 2 . (80) 

Using the steady-state solution for vL (equation |70j>, we equate Q + to the rate of 
local energy loss by radiation. If the disk is optically thick, the disk radiates (from 
both sides) at a rate Q = 2cr7’jj sk , with cr being the Stefan-Boltzmann constant. 
This yields an effective temperature profile. 


4 _3GM,M 
disk 8kg 



(81) 


Away from the inner boundary, the steady-state temperature profile for a viscous 
disk (7dj s k oc r~ 3 / 4 ) is steeper than for irradiation. For any accretion rate, we then 
expect viscous heating to be most important in the inner disk, whereas irradiation 
always wins out at sufficiently large radii. 

The viscous disk temperature profile is not what we get from considering just the 
local dissipation of potential energy. The gradient of the potential energy per unit 
mass e, is de/d r = GM^/r 1 . For an accretion rate M, the luminosity available to be 
radiated from an annulus of width Ar due to local potential energy release would 
be. 


L = 


1 GM,M 

2 r 2 


Ar, 


(82) 
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where the factor of a half accounts for the fact that half the energy goes into in¬ 
creased kinetic energy, with only the remainder available to be thermalized and 
radiated. Equating this luminosity to the black body emission from the annulus, 
2°'7(i!sk ‘ 2 nrAr, would give a profile that is a factor three different from the asymp¬ 
totic form of equation ( |8T| . The difference arises because the radial transport of 
angular momentum is accompanied by a radial transport of energy. The local lu¬ 
minosity from the disk surface at any radius then has a contribution from potential 
energy liberated closer in. 

The optically thick regions of irradiated protoplanetary disks will be vertically 
isothermal. When viscous heating dominates, however, there must be a vertical tem¬ 
perature gradient to allow energy to be transported from the mid-plane toward the 
photosphere. What this gradient looks like, in detail, depends on the vertical distri¬ 
bution of the heating, which is not well known. However, an approximation to T (z) 
can be derived assuming that the energy dissipation due to viscosity is strongly con¬ 
centrated toward the mid-plane. We define the optical depth to the disk mid-plane, 

t=^k r E, (83) 


where Kr is the Rosseland mean opacity. The vertical density profile of the disk is 
p(z). If the vertical energy transport occurs via radiative diffusion then for T 1 
the vertical energy flux F(z) is given by the equation of radiative diffusion 12991 


F z {z) 


16ar 3 dr 
3 K R p dz ’ 


(84) 


Now assume for simplicity that all of the dissipation occurs at z = 0. In that case 
F z (z) = crTjj sk is constant with height. We integrate from the mid-plane to the pho¬ 
tosphere at z p h assuming that the opacity is also constant. 
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(85) 

( 86 ) 


where the final equality relies on the fact that for T> 1 almost all of the disk gas 
lies below the photosphere. For large optical depth 7' L 4 7’jj sk and the equation 
simplifies to. 



(87) 


Often both stellar irradiation and accretional heating contribute significantly to the 
thermal balance of the disk. If we define 7di s k,visc to be the effective temperature 
that would result from viscous heating in the absence of irradiation (i.e. the quantity 
called T(ji s k, with no subscript, above) and T m to be the irradiation-only effective 
temperature, then. 
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^^<v, c + ^ r (88) 

is an approximation for the central temperature, again valid for T 3g> 1. 

These formulae can be applied to estimate the location of the snow line. In the So¬ 
lar System, meteoritic evidence 112381 places the transition between water vapor and 
water ice, which occurs at a mid-plane temperature of 150-180 K, at around 2.7 AU. 
This is substantially further from the Sun than would be expected if the only disk 
heating source was starlight. Including viscous heating, however, an accretion rate 
of 2 x 10 8 Mq yr 1 could sustain a mid-plane temperature of 170 K at 2.7 AU in a 
disk with E = 400 g cm -2 and Kr= 1 cm 2 g 1 . This estimate (from equation|87|) is 
consistent with more detailed models for protoplanetary disks (491 . though consid¬ 
erable variation in the location of the snow line is introduced by uncertainties in the 
vertical structure (227). 


3.4 Warped disks 


The classical equation for surface density evolution needs to be rethought if the disk 
is non-planar. Disks may be warped for several reasons; the direction of the angular 
momentum vector of the gas that forms the disk may not be constant, the disk may 
be perturbed tidally by a companion 1 192 2461 . or warped due to interaction with 
the stellar magnetosphere 11911 , 

A warp affects disk evolution through physics that is independent of its origin 
(for a brief review, see 12451 ). In a warped disk, neighboring annuli have specific 
angular momenta that differ in direction as well as in magnitude. If we define a unit 
tilt vector \{r.t) that is locally normal to the disk plane, the shear then has a vertical 
as well as a radial component (2491 , 


S = 


d n <91 

r-l + rQ^-. 
dr dr 


(89) 


The most important consequence of the vertical shear is that it introduces a periodic 
vertical displacement of radially separated fluid elements. As illustrated in Figure[TT] 
this displacement, in turn, results in a horizontal pressure gradient that changes sign 
across the mid-plane and is periodic on the orbital frequency. In a Keplerian disk 
this forcing frequency is resonant with the epicyclic frequency. 

How the disk responds to the warp-generated horizontal forcing depends on the 
strength of dissipation (263) . If the disk is sufficiently viscous, specifically if, 

a > (90) 

r 

the additional shear is damped locally. The equation for the surface density and tilt 
evolution (the key aspects of which are derived in (278) . though see (247) for a 
complete treatment) then includes terms which diffusively damp the warp at a rate 
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d p|dv 


Fig. 11 Illustration (after ||2 08ll2451 ) of how a warp introduces an oscillating radial pressure gra¬ 
dient within the disk. As fluid orbits in a warped disk, vertical shear displaces the mid-planes of 
neighboring annuli. This leads to a time-dependent radial pressure gradient dP/dr(z). Much of the 
physics of warped disks is determined by how the disk responds to this warp-induced forcing. 


that is related to the radial redistribution of angular momentum. Normally, warp 
damping is substantially faster than the viscous evolution of a planar disk. Even 
for a Navier-Stokes viscosity — which is fundamentally isotropic — the effective 
viscosity which damps the warp is a factor sa 1 /2a 2 larger than its equivalent in a 
flat disk. Rapid evolution also occurs in the opposite limit of an almost inviscid disk 
with, 

(x<~, (91) 

r 

but in this case the component of angular momentum associated with the warp is 
communicated radially in the form of a wave. For a strictly inviscid disk the lin¬ 
earized fluid equations for the evolution of the tilt vector take a simple form 12151 . 


dh 

dt 2 


i d 

Er 2 Q dr 




(92) 


The speed of the warp wave is v w ~c s / 2. 

In most cases we expect protoplanetary disks to have a <h/r, and warps will 
evolve in the wave-like regime. We expect, however, that the details of warp 
evolution will depend upon the nature of angular momentum transport, and 
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Fig. 12 Illustration, after Spruit ( 3138 . of the different regions of a disk wind solution. 


nothing is known about how warps behave for the transport mechanisms (such 
as non-ideal MHD) most relevant to protoplanetary disks. The differences be¬ 
tween warp evolution with realistic transport, and that with a Navier-Stokes 
viscosity, are of undetermined size. 


3.5 Disk winds 

Viscous evolution driven by redistribution of angular momentum is a consequence 
of turbulent (or possibly laminar) stresses that are internal to the fluid. Evolution 
can also be driven by external torques, the most important of which is the magnetic 
torque that an MHD wind exerts on the surface of the disk. An excellent pedagogical 
introduction to disk winds is the review by Spruit 1131311 . while Konigl & Salmeron 
11801 provide a more recent account that addresses the peculiarities specific to pro¬ 
toplanetary disks. 

Winds that are driven solely by pressure gradients (“thermal winds”) are of in¬ 
terest as a mechanism for disk dispersal ((J8| but do not otherwise change the qual¬ 
itative character of disk evolution. Instead, we consider here MHD winds, with the 
simplest case being a well-ionized disk is threaded by a large-scale ordered poloidal 
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magnetic field lf57l . In the ideal MHD limit the fluid is tied to magnetic field lines, 
which can facilitate acceleration of the wind while exerting a back reaction on the 
disk that removes angular momentum. 

Blandford-Payne winds are not the only type of MHD outflow of interest for 
protoplanetary disk systems. Outflows could be launched by a gradient of 
toroidal magnetic field pressure |21X|. perhaps during or shortly after the col¬ 
lapse of the cloud that forms both the star and the disk. Jets can also originate 
from the interaction between the stellar magnetosphere and the disk 113051 . 

The structure of such a wind, illustrated in Figure [12J generically has three re¬ 
gions. Within the disk the energy density in the magnetic field, B 2 / 8k, is smaller 
than pc 2 , the thermal energ\fj Due to flux conservation, however, the energy in the 
vertical field component, Bz /8k, is roughly constant with height for z < r, while 
the gas pressure typically decreases at least exponentially with a scale height K<r 
This leads to a region above the disk surface where magnetic forces dominate. The 
magnetic force per unit volume is. 


where the current, 


JxB 


B 2 \ B VB 

87T/ + 4k 

(93) 

- V x B. 

(94) 


The force, in general, can be written as shown above as the sum of a magnetic 
pressure gradient and a force due to magnetic tension. In the disk wind region where 
magnetic forces dominate, the requirement that they exert a finite acceleration on the 
low density gas implies that the force approximately vanishes, i.e. that. 


JxBssO. 


(95) 


The stmcture of the magnetic field in the magnetically dominated region is then 
described as being “force-free”, and in the disk wind case (where B changes slowly 
with z) the field lines must be approximately straight to ensure that the magnetic 
tension term is also individually small. If the field lines support a wind, the force- 
free structure persists up to where the kinetic energy density in the wind, pv 2 , first 
exceed the magnetic energy density. This is called the Alfven surface. Beyond the 
Alfven surface, the inertia of the gas in the wind is sufficient to bend the field lines, 
which tend to wrap up into a spiral structure as the disk below them rotates. 

Magneto-centrifugal driving can launch a wind from the surface of a cold gas 
disk if the magnetic field are sufficiently inclined to the disk normal. The critical 
inclination angle in ideal MHD can be derived via an exact mechanical analogy. To 


3 We can also consider situations where the magnetic pressure in the disk is stronger than the gas 
pressure, though it must always be weaker than pv^. 
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Fig. 13 Geometry for the calculation of the critical angle for magneto-centrifugal wind launching. 
A magnetic field line s, inclined at angle 6 from the disk normal, enforces rigid rotation at the 
angular velocity of the foot point, at cylindrical radius tx! = ro in the disk. Working in the rotating 
frame we consider the balance between centrifugal force and gravity. 


proceed, we note that in the force-free region the magnetic field lines are (i) basically 
straight lines, and (ii) enforce rigid rotation out to the Alfven surface at an angular 
velocity equal to that of the disk at the field line’s footpoint. The geometry is shown 
in Figure [13] We consider a field line that intersects the disk at radius ro, where the 

angular velocity is f2o = GM , : /and that makes an angle 0 to the disk normal. 
We define the spherical polar radius r , the cylindrical polar radius QJ, and measure 
the distance along the field line from its intersection with the disk at z = 0 as s. In 
the frame co-rotating with Qq there are no magnetic forces along the field line to 
affect the acceleration of a wind; the sole role of the magnetic field is to constrain 
the gas to move along a straight line at constant angular velocity. Following this 
line of argument, the acceleration of a wind can be fully described in terms of an 
effective potential, 

‘S’effO) = _ If2 ( }G7 2 (.s). (96) 

r(s) 2 

The first term is the gravitational potential, while the second describes the centrifu¬ 
gal potential in the rotating frame. 

Written out explicitly, the effective potential is. 


‘J’efftA) 


_ GM *_ 

(s 2 + 2sro sin 0 + Tq) */ 2 


1 

2 


Q.q (ro + ssin0) 2 . 


(97) 


This function is plotted in Figure[l4]for various values of the angle 0. For a vertical 
field line (0=0) the effective potential is a monotonic ally increasing function of 
distance ,v, for modest values of 0 there is a potential barrier defined by a maximum 
at some s = .v ma x. while for large enough 0 the potential decreases from ,y = 0. In this 
last case purely magneto-centrifugal forces suffice to accelerate a wind off the disk 








38 


Philip J. Armitage 



Fig. 14 The variation of the disk wind effective potential 4> e ff (in arbitrary units) with distance s 
along a field line. From top downwards, the curves show field lines inclined at 0°, 10°, 20°, 30° (in 
bold) and 40° from the normal to the disk surface. For angles of 30° and more from the vertical, 
there is no potential barrier to launching a cold MFID wind directly from the disk surface. 


surface, even in the absence of any thermal effects. The critical inclination angle of 
the field can be found by computing 0 cr j t , specified though the condition. 


<9 2 <J> 


'eff 


ds 2 


= 0 . 


s—0 


Evaluating this condition, we find, 

1-4 sin 2 0 cr j t = 0 
=► flcrit = 30°, 


(98) 


(99) 


as the minimum inclination angle from the vertical needed for unimpeded wind 
launching in ideal MHD (57). Since most of us are more familiar with mechanical 
rather than magnetic forces, this derivation in the rotating frame offers the easiest 
route to this result. But it can, of course, be derived just as well by working in the 
inertial frame of reference (3131 . 

The rigid rotation of the field lines interior to the Alfven surface means that gas 
being accelerated along them increases its specific angular momentum. The mag¬ 
netic field, in turn, applies a torque to the disk that removes a corresponding amount 
of angular momentum. If a field line, anchored to the disk at radius ro, crosses the 
Alfven surface at (cylindrical) radius r\, it follows that the angular momentum flux 
is. 
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L w — . 


( 100 ) 


where M w is the mass loss rate in the wind. Removing angular momentum at this 
rate from the disk results in a local accretion rate M = 2L w /Qo r o■ The ratio of the 
disk accretion rate to the wind loss rate is, 



( 101 ) 


If r,\ substantially exceeds rg (by a factor of a few, which is reasonable for detailed 
disk wind solutions) a relatively weak wind can carry away enough angular momen¬ 
tum to support a much larger accretion rate. 

The behavior of a disk that evolves under wind angular momentum loss depends 
on how the wind and the poloidal magnetic field respond to the induced accretion. 
It is not immediately obvious that a steady accretion flow is even possible. The 
form of the effective potential (Figure [l4| suggests that the rate of mass and angular 
momentum loss in the wind ought to be a strong function of the inclination of the 
field lines — for 0 < 30° there is a potential barrier to wind launching, while for 
6 > 30° there is no barrier at all. How 6 responds to changes in the inflow rate 
through the disk is of critical importance 1217 66] 1250 1. and there is no simple 
analog of the diffusive disk evolution equation. 

Irrespective of the (uncertain) details, viscous and wind-driven disks exhibit 
some qualitative difference that may enable observational tests. The classical 
test is the evolution of the outer disk radius, which expands in viscous models 
(if there is no mass loss, even in the form of a thermal wind) but contracts 
if an MHD wind dominates. Old and almost forgotten observations f3T21 of 
disk radius changes in dwarf novae (accreting white dwarfs in mass transfer 
binary systems) provided empirical support for viscous disk evolution in those 
specific systems. Disk winds also remove energy, and so another potential test 
is to look for evidence of the dissipation of accretion energy within the disk 
that is present in viscous models but absent for winds. 


3.5.1 Magnetic field transport 

The strength and radial profile of the vertical magnetic field threading the disk are 
important quantities for disk winds, and for turbulence driven by MHD processes. 
Disks form from the collapse of magnetized molecular clouds, and it is inevitable 
that they will inherit non-zero flux at the time of formation. The poloidal component 
of that flux can subsequently be advected radially with the disk gas, diffuse relative 
to the gas, or (if the flux has varying sign across the disk) reconnect. 

A theory for the radial transport of poloidal flux within geometrically thin ac¬ 
cretion disks was developed by Lubow, Papaloizou & Pringle (1994) 12161 . They 
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considered a disk within which turbulence generates an effective viscosity v and an 
effective magnetic diffusivity 77 . The disk is threaded by a vertical magnetic field 
B z (r,t), which is supported by a combination of currents within the disk and (po¬ 
tentially) a current external to the disk. Above the disk, as in Figure [12] the field is 
force-free. The field lines bend within the disk, such that the poloidal field has a 
radial component B rs (r,t) at the disk surface. 

To proceed (following the notation in 1 134 1), we first write the poloidal field in 
terms of a magnetic flux function t//, such that B = V i/r x e,.,, where is a unit 
vector in the azimuthal directions. The components of the field are. 


Br 

B z 


1 dy 
r dz ’ 
1 dy 
r dr 


( 102 ) 


From the second of these relations we note that y is (up to a factor of 2n) just 
the vertical magnetic flux interior to radius r. We split y into two pieces, a piece 
t/Zdisk due to currents within the disk, and a piece y*, due to external currents (“at 
infinity”), 

W= V / disk + V / ~- (103) 

The external current generates a magnetic field that is uniform across the disk. 

With these definitions, the evolution of the poloidal field in the simplest analysis 
C36l[l34l obeys, 

^ + r (v adv B z + v diff B„) = 0, (104) 

where v adv is the advective velocity of magnetic flux and i’,iifr its diffusive velocity 
due to the turbulent resistivity within the disk. The disk component of the flux func¬ 
tion is related to the surface radial field via an integral over the disk. Schematically, 

ft 'out 

VAiisk(r) = / F{rJ)B rs {r')&r', (105) 

where F is a rather complex function that can be found in Guilet & Ogilvie (2014) 
1134!]. The appearance of this integral reflects the inherently global nature of the 
problem — a current at some radius within the disk affects the poloidal magnetic 
field everywhere, not just locally — and makes analytic or numerical solutions for 
flux evolution more difficult. Nonetheless, equation ( |105| > can be inverted to find B rs 
from y, after which the more familiar equation ( |104| ) can be solved for specified 
transport velocities to determine the flux evolution. 

Equation ( |104| > expresses a simple competition, the inflow of gas toward the star 
will tend to drag in poloidal magnetic field, but this will set up a radial gradient 
and be opposed by diffusion. The physical insight of Lubow et al. (1994) 12161 was 
to note that although both of these are processes involving turbulence (and, very 
roughly, we might guess that v ~ !]), the scales are quite distinct. From Figure [12] 
we note that because the field lines bend within the disk, a moderately inclined 
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external field (with B rs ~ B-) above the disk only has to diffuse across a scale ~ h to 
reconnect with its oppositely directed counterpart below the disk. Dragging in the 
field with the mean disk flow, however, requires angular momentum transport across 
a larger scale r. In terms of transport velocities, in a steady-state we have, 


' ; adv 


t'diff 


V 
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r 

V Brs 
h B z ’ 


(106) 


so that for v ~ q and B rs ~ B- diffusion beats advection by a factor ~ ( h/r)~ l 
1. Defining the magnetic Prandtl number P m = v/rj as the ratio of the turbulent 
viscosity to the turbulent resistivity, we would then expect that in steady-state 12161 . 


B n 

B- 



(107) 


This argument is the origin of the claim that thin disks do not drag in external 
magnetic fields. It suggests that obtaining enough field line bending to launch a 
magneto-centrifugal wind ought to be hard, and that whatever primordial flux the 
disk is born with may be able to escape easily. 

The physical arguments given above are robust, but a number of authors have em¬ 
phasized that the calculation of the transport velocities that enter into equation (104 1 
involves some subtleties 112501 [2 13 ;, [ 133M 34 1323 1. The key point is that the viscos¬ 
ity and resistivity that enter into the equation for flux transport should not be com¬ 
puted as density-weighted vertical averages, but rather (in the case of the induction 
equation) as conductivity-weighted averages 11331 . This makes a large difference 
for protoplanetary disks, where the conductivity is both generally low, and highest 
near the disk surface where the density is small. The derived transport velocities 
are, moreover, functions of the poloidal field strength, in the sense that diffusion 
becomes relatively less efficient as the field strength decreases. It should be noted 
that none of the flux transport calculations fully includes all of the MHD effects ex¬ 
pected to be present in protoplanetary disks (see (4.4.3 i. It seems possible, though, 
that the variable efficiency of flux diffusion could simultaneously allow. 


• For rapid flux loss from the relatively strongly magnetized disks formed from star 
formation 12031 . averting overly rapid wind angular momentum loss that would 
be inconsistent with observed disk lifetimes. 

• For convergence toward a weak but non-zero net poloidal flux (possibly with a 
ratio of thermal to poloidal field magnetic pressure at the mid-plane ji ~ 10 4 — 
10 7 ) later in the disk lifetime Ml 341 . 


As we will discuss in the next section, poloidal field strengths in roughly this range 
are of interest for their role in stimulating MHD instabilities within weakly ionized 
disks, so this is a speculative but interesting scenario. 
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4 Turbulence 


Turbulence within protoplanetary disks is important for two independent reasons. 
First, if it is strong enough and has the right properties, it could account for disk 
evolution by redistributing angular momentum much faster than molecular viscos¬ 
ity. Second, turbulence has its fingers in a plethora of planet formation processes, 
ranging from the collision velocities of small particles 12531 to the formation of 
planetesimals lll60ft and the migration rate of low-mass planets ll77l . For these rea¬ 
sons we would like to understand disk turbulence, even if (as is possible) it is not 
always responsible for disk evolution. 

The first order of business when considering possibly turbulent fluid systems is 
usually to estimate the Reynolds number, which is a dimensionless measure of the 
relative importance of inertial and viscous forces. For a system with characteristic 
size L, velocity U , and (molecular) viscosity v, the Reynolds number is defined. 


Re = 


UL 


(108) 


There is no unique or “best” definition of U and L for protoplanetary disks, but 
whatever choice we make gives a very large number. For example, taking L = li and 
JJ = c s then our estimate of the viscosity at 1 AU (equation 62 1 implies Re ~ 10 11 . 


By terrestrial standards this is an enormous Reynolds number. Experiments on flow 
through pipes, for example — including those of Osborne Reynolds himself — show 
that turbulence is invariably present once Re > 10 4 (98). If turbulence is present 
within protoplanetary disks, there is no doubt that viscous forces will be negligible 
on large scales, and the turbulence will exhibit a broad inertial range. 

Figure [l5]lists some of the possible sources of turbulence within protoplanetary 
disks. It’s a long list! We can categorize the candidates according to various criteria, 

• The physics involved in generating the turbulence. The simplest possibility 
(which appears unlikely) is that turbulence develops spontaneously in an isother¬ 
mal, purely hydrodynamical shear flow. More complete physical models invoke 
entropy gradients, disk self-gravity or magnetic fields as necessary elements for 
the origin of turbulence. 

• The origin of the free energy that sustains the turbulence, which could be the 
radial or vertical shear, heating from the star, or velocity differences between gas 
and solid particles. 

• The character of the instabilities proposed to initiate turbulence from an ini¬ 
tially non-turbulent flow. Linear instabilities grow exponentially from arbitrar¬ 
ily small perturbations, while non-linear instabilities require a finite amplitude 
disturbance. Demonstrating the existence of linear instabilities is relatively easy, 
whereas proving that a fluid system is non-linearly stable is very hard. 

• The species involved. In this section we concentrate on instabilities present in 
purely gaseous disks; additional instabilities are present once we consider how 


gas interacts aerodynamic ally with its embedded solid component (16.3 i. 
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Figure [16] illustrates the dominant fluid motions or forces involved in some of the 
most important disk instabilities. 

For each candidate instability we would like to know the disk conditions under 
which it would be present, its growth rate, and the strength and nature of the tur¬ 
bulence that eventually develops. For disk evolution we are particularly interested 
in how efficiently the turbulence transports angular momentum (normally charac¬ 
terized by an effective a). In most cases the efficiency of transport can only be 
determined using numerical simulations, whose fluctuating velocity and magnetic 
fields can be analyzed to determine a via the relation l38l . 


a = 


SvrSvij) 


c 


2 

s 


B r B(j) \ 

^P c2 S Ip ’ 


(109) 


where the angle brackets denote a density weighted average over space (and time, 
in some instances). The first term in this expression is the Reynolds stress from 
correlated fluctuations in the radial and perturbed azimuthal velocity, the second 
term in the Maxwell stress from MHD turbulence. We speak of the stress as being 
“turbulent” if the averages in the above relation are dominated by contributions from 
small spatial scales. It is also possible for a disk to sustain large scale stresses — 
for example at some radius we might have non-zero mean radial and azimuthal 
magnetic fields — which are normally described as being “laminar”. 
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Fig. 16 A summary of the most important instabilities that can be present in protoplanetary disks. 
Self-gravity is important for sufficiently massive and cold disks. It leads to spiral arms and gravita¬ 
tional torques between regions of over-density. The magnetorotational instability occurs whenever 
a weak magnetic field is sufficiently coupled to differential rotation. The magnetic field acts to cou¬ 
ple fluid elements at different radii, leading to an instability that can sustain MHD turbulence and 
angular momentum transport. The vertical shear instability feeds off the vertical shear that is set 
up in disks with realistic temperature profiles. It is a linear instability characterized by near-vertical 
growing modes. The subcritical baroclinic instability is a non-linear instability that operates in the 
presence of a sufficiently steep radial entropy gradient. It resembles radial convection, and leads to 
self-sustained vortices within the disk. 


4.1 Hydrodynamic turbulence 


The dominant motion in protoplanetary disks is Keplerian orbital motion about a 
central point mass. Simplifying as much as possible, we first ask whether, in the 
absence of magnetic fieldtjj the radial shear present in a low-mass disk would be 
unstable to the development of turbulence. We first consider (rather unrealistically) 


a radially isothermal disk, where according to equation (17 1 there is no vertical 
shear. We then turn to the more general case where the temperature varies with 
radius, giving rise both to vertical shear and qualitatively distinct possibilities for 
instability. 


4 Ignoring magnetic fields in astrophysical accretion flows is generally a stupid thing to do, and 
indeed there is broad consensus that the magnetorotational instability (MRI) CD is responsible 
for turbulence and angular momentum transport in most accretion disks. In protoplanetary disks, 
however, the low ionization fraction means that the dominance of MHD instabilities is much less 
obvious, and purely hydrodynamic effects could in principle be important. 
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4.1.1 Linear and non-linear stability 

The linear stability of a shear flow with a smoothly varying Q (r) against axisym- 
metric perturbations is given by Rayleigh’s criterion (this is derived in most fluids 
textbooks, see e.g. 12791 ’). The flow is stable if the specific angular momentum in¬ 
creases with radius, 

— = — (r 2 £ 2) >0-4 stability. (110) 

dr dr v ' 

A Keplerian disk has l ^Jr and is linearly stable. 

There is no mathematical proof of the non-linear stability of Keplerian shear 
flow, but nor is there any known instability. The apparently analogous cases of pipe 
flow and Cartesian shear flows — which are linearly stable but undergo non-linear 
transitions to turbulence — are in fact sufficiently different problems as to offer no 
guidance ll38l . There are analytic and numerical arguments against the existence 
of non-linear instabilities li39l , which although not decisive [ 2871 essentially rule 
out the hypothesis that a non-linear instability could result in astrophysically in¬ 
teresting levels of turbulence 11961 . The same conclusion follows from laboratory 
experiments that have studied the stability of quasi-Keplerian rotation profiles in 
Taylor-Couette experiments |[99l . 

Laboratory experiments, and most theoretical work, consider the stability of 
unstratified cylindrical shear flows. It has been suggested that new instabili¬ 
ties (of a distinct character, related to the existence of locations in the flow 
known as critical layers, for a review see 12311 ) arise when the vertical strati¬ 
fication present in disks is included B2251 . Study of this possibility remains in 
its infancy. 


4.1.2 Entropy-driven instabilities 

A separate class of purely hydrodynamic instabilities (no self-gravity, no mag¬ 
netic fields) are what might loosely be called “entropy-driven” instabilities, in that 
they rely on the existence of a non-trivial temperature structure. The prototypical 
entropy-driven instability is of course convection, which could occur in the vertical 
direction if dissipation (associated with the physical process behind angular mo¬ 
mentum transport) sets up an unstable entropy profile. This is evidently only con¬ 
ceivable in the region where viscous dissipation dominates, as irradiation prefers 
a nearly isothermal vertical structure. Even there, convective turbulence in disks is 
less efficient at transporting angular momentum than it is in transporting heat 1119711 . 
and this disparity creates a formidable barrier to creating consistent models in which 
convection is the primary source of disk turbulence. Convection may still be present 
in some regions of disks, perhaps especially at high accretion rates, but as a byprod- 
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uct of independent angular momentum transport processes (for an example in dwarf 
novae, see GH). 

A disk that has a radial temperature gradient necessarily has vertical shear (equa¬ 
tion 17 1 . The free energy associated with the vertical shear can be accessed via 


the vertical shear instability (VSI) analyzed by Nelson, Gressel & Umurhan (2013) 
1243 1 . The VSI is a disk application of the Goldreich-Schubert-Fricke instability 
1 1261 . 1_12] of rotating stars, and was proposed as a source of protoplanetary disk 
transport by Urpin & Brandenburg 038 The VSI is a linear instability with a 
maximum growth rate that is of the order of IiQk 112061 . but which is strongly de¬ 
pendent on the radiative properties of the disk. The reason is that to access the free 
energy in the vertical shear requires vertical fluid displacements, which are easy in 
the limit that the disk is strictly vertically isothermal but strongly suppressed if it is 
stably stratified. The local cooling time of the fluid is thus a critical parameter, and 
the VSI will only operate in regions of the disk where radiative cooling and heating 
processes result in a cooling time that is the same or shorter than the dynamical time 
.Qjf 1 . This, in practice, restricts the application of the VSI to intermediate radii (Lin 
& Youdin suggest 5-50 AU 120611 1. and limits its effectiveness if the dust opacity is 
reduced (due to coagulation into large particles). Under the right conditions, how¬ 
ever, numerical simulations suggest that the VSI can generate relatively small but 
possibly significant levels of transport, with both Nelson et al. 12431 and Stoll & 
Kley 13181 finding a of a few x 10 -4 . 

The radial entropy gradient may itself be unstable. The simplest instability would 
be radial convection (a linear instability). For a disk with pressure profile P(r), den¬ 
sity profile p (r), and adiabatic index y, we define the Brunt-Vaisala frequency. 


Nt 


1 dP d / P \ 
yp dr dr \p r ) 


( 111 ) 


The Solberg-Hoi'land criterion indicates that a Keplerian disk is convectively unsta¬ 
ble if, 

N? + nl< 0 . ( 112 ) 

Protoplanetary disks never (or at least almost never) have a steep enough profile of 
entropy to meet this condition, so radial convection will not set in. A different insta¬ 
bility (the subcritical baroclinic instability, SBI) is possible, however, if the weaker 
condition V, 2 < 0 (which is just the Schwarzschild condition for non-rotating con¬ 
vection) is satisfied 1 2681,1269 4991 1. The SBI, which is likely related to observations 
of vortex formation in earlier numerical simulations 111751 . is a non-linear instability 
that can be excited by finite amplitude perturbations. (Confusingly, it is unrelated to 
the linear “baroclinic instability” studied in planetary atmospheres.) The SBI relies 
on radial buoyancy forces to sustain vortical motion via baroclinic driving. This type 
of effect is possible in disks in which surfaces of constant density are not parallel to 
surface of constant pressure. Mathematically, for a fluid with vorticity ft) = V x v, 
we can take the curl of the momentum equation to get an equation for the vortensity 


5 As we shall see, a general rule is that all disk instabilities have long histories and pre-histories. 
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co/p. 


D f co 


Dr \p 



(113) 


The baroclinic term, which for the SBI is responsible for generating and maintain¬ 
ing vorticity in the presence of dissipation, is the second term on the right hand side. 
The SBI, as with the VSI, is sensitive to the cooling time Ml99112831 , in this case 
because the baroclinic driving depends on the disk neither cooling too fast (which 
would eliminate the buoyancy effect) nor too slow (which would lead to constant 
temperature around the vortex). In compressible simulations, Lesur & Papaloizou 
(2010) Il99tt found that under favorable disk conditions the SBI could lead to out¬ 
ward transport of angular momentum with a ~ 10 1 . 


4.2 Self-gravity 


A disk is described as self-gravitating if it is unstable to the growth of surface den¬ 
sity perturbations when the gravitational force between different fluid elements in 
the disk is included along with the force from the central star. For a disk with sound 
speed c s , surface density E and angular velocity £2 (assumed to be close to Keple- 
rian) a linear analysis (for textbook treatments, see e.g. II141 12791 ) shows that a disk 
becomes self-gravitating when the Toomre Q 13281 . 


„ c s £2 „ 

Q = < Cork, 

nLrL 


(114) 


where (9 cnt ~ 1. We can deduce this result informally using an extension of the time 
scale argument that gives the thermal Jeans mass. We first note that pressure will pre¬ 
vent the gravitational collapse of a clump of gas, on scale Ar. if the sound-crossing 
time Ar/c s is shorter than the free-fall time \JAr 2 /GAr 2 E. (We’re ignoring factors 
of 2, 71 and so on.) Equating these time scales gives the minimum scale that might 
be vulnerable to collapse as Ar ~ c 2 s /GE. On larger scales, collapse can be averted 
if the free-fall time is longer than the time scale on which radial shear will separate 
initially neighboring fluid elements. For a Keplerian disk, this time scale is ~ O' 1 . 
If the disk is just on the edge of instability, the minimum collapse scale set by pres¬ 
sure support must equal the maximum collapse scale set by shear. Imposing this 
condition for marginal stability we obtain c s £2/GE ~ 1, in accord with the formal 
result quoted above. 

To glean some qualitative insight into where a disk might be self-gravitating, 
consider a steady-state disk that is described by an a model in which the transport 
arises from some process other than self-gravity. Collecting some previous results, 
the steady-state condition implies vE = M/3ti, the a prescription is v = ac s h. and 


hydrostatic equilibrium gives h = c s /£ 2. Substituting into equation (114) we find, 


Q = 


3 ac^ 
GM ' 


(115) 
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Protoplanetary disks generically get colder (and hence have lower c s ) at larger dis¬ 
tances from the star, and this is where self-gravity is most likely to be important. 

The disk mass required for self-gravity to become important can be estimated. 
Ignoring radial gradients of all quantities, we write the disk mass M^sk ~ 7ii~E, and 
again use the hydrostatic equilibrium result Ir = c s /£ 2. Equation (114 1 then gives. 


^disk 

M* 



(116) 


as the condition for instability. This manipulation of a local stability criterion into 
some sort of global condition is ugly, and begs the question of where in the disk 
Mdisk and h/r should be evaluated. We can safely conclude, nonetheless, that for a 
typical protoplanetary disk with (h/r) ~ 0.05 a disk mass of 10~ 2 M* will not be 
self-gravitating, whereas one with 0.1 M* may well be. 

There are two possibles outcomes of self-gravity in a disk, 

• The disk may establish a (quasi) stable state, characterized globally by trailing 
spiral overdensities. Gravitational torques between different annuli in the disk 
transport angular momentum outward, leading to accretion. 

• The pressure and tidal forces, which by definition are unable to prevent the onset 
of gravitational collapse, may never be able to stop it once it starts. In this case 
the disk fragments into bound objects, which interact with (and possibly accrete) 
the remaining gas. 

Both possibilities are of interest. Angular momentum transport due to self-gravity 
may be dominant, at least on large scales, at early times while the disk is still mas¬ 
sive. Fragmentation, which was once considered a plausible mechanism for form¬ 
ing the Solar System’s giant planets DM), remains of interest as a way to form 
sub-stellar objects and (perhaps) very massive planets. We will look at both in turn. 

Gravity is a long-range force, and it is not at all obvious that we can deploy the 
machinery developed for viscous disks to study angular momentum transport in a 
self-gravitating disk. The transport could be largely non-local, driven for example 
by large-scale structures in the density field (such as bars) or by waves that transport 
energy and angular momentum a significant distance before dissipating |40). There 
is no precise criterion for when self-gravitating transport can be described using a 
local theory, but numerical simulations indicate that this is a reasonable approxima¬ 
tion for low-mass disks with « 0.1 1 209,,86 . 109,1314 1. Transport in more 

massive disks, such as might be present during the Class 0 and Class I phases of star 
formation, cannot be described locally (for multiple reasons, e.g. 1210 13311 ). 

Assuming that a local description of the transport is valid, we can use a ther¬ 
mal balance argument to relate the efficiency of angular momentum transport to the 
cooling time. Adopting a one-zone model for the vertical structure, we define the 
thermal energy of the disk, per unit surface area, as, 
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where c s is the mid-plane sound speed and 7 is the adiabatic index. The cooling time 
(analogous to the Kelvin-Helmholtz time for a star) is then, 


'cool — 


u 


2a T 4 ’ 

ZOi disk 


(118) 


where 7 disk is the effective temperature. Equating the cooling rate, Q . = 2oT d 4 to 
the local viscous heating rate, Q+ = (9/4 )vEQ 2 (equation |80j), and adopting the 
a-prescription (equation [65]), we find. 


4 1 

97(7-1) Qt coo i ' 


(119) 


This relation, which is a general property of a disks quite independent of self¬ 
gravity, just says that a rapidly cooling disk needs efficient angular momentum 
transport if it to generate heat fast enough to remain in thermal equilibrium. 

For most sources of angular momentum transport we are no more able to de¬ 
termine f coo i from first principles than we are a, so the above relation does not 
move us forward. Self-gravitating disks, however, have the unusual property that 
their Toomre Q, measured in the saturated (non-linear) state, is roughly constant 
and similar to the critical value Qcrit determined from linear theory. This property 
arises, roughly speaking, because the direct dependence of the linear stability crite¬ 
rion on temperature (via c s ) invites a stabilizing feedback loop — a disk that cools 
so that Q < Qcnt is more strongly self-gravitating, and produces more heating, while 
one that heats so that Q > (/ tnt shuts off the instability. It is therefore reasonable to 
assume that a self-gravitating disk that does not fragment maintains itself close to 
marginal stability, as conjectured by Paczynski (1978) 12611 . 

If we assume that Q = Qo exactly (where Q 0 is some constant presumably close 
to 2crit) then we have enough constraints to explicitly determine the functional form 
of a for a self-gravitating disk. Since Q depends on the mid-plane sound speed, 
c s = sJksTc/ the condition of marginal stability directly gives us T C (L.Q), 

r t = *w(^)g. (120) 


We can use this to determine f co 0 i, and from that a, with the aid of the vertical 
structure relations developed in (3.3 To keep things simple, we adopt an opacity. 


k r = k 0 T, 2 


( 121 ) 


that is appropriate for ice grains, and assume the disk is optically thick. The opacity 
law, together with the relations for the optical depth, T = (1/2 )Ekr, and the mid¬ 
plane temperature, 7/ 4 /7’ d i j sk ~ (3/4 )t, then leads to, 


a 


_64n 2 QlG 2 o 

27 Kq \k B ) 


2 


( 122 ) 
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which coincidentally (for this opacity law) is only a function of Q. It may look 
cumbersome — and the numerical factors are certainly not to be trusted — but what 
we have shown is that for a locally self-gravitating disk a is simply a constant times 
a determined function of E and Q. This result allows for the evolution of low-mass 
self-gravitating disks to be modeled as a pseudo-viscous process 112001 [781 [284 1. 

Self-gravity is typically important in protoplanetary disks at large radii, where ir¬ 
radiation is usually the dominant factor determining the disk’s thermal state (except 
at high accretion rates). The generalization of the self-regulation argument given 
above is obvious; if irradiation is not so strong as to stabilize the disk on its own 
then viscous heating from the self-gravitating “turbulence” has to make up the dif¬ 
ference. The partially irradiated regime of self-gravitating disks has been studied 
using local numerical simulations 112891 . and the analytic generalization for the ef¬ 
fective a that results can be found in Rafikov (2015) 12851 . 

Ignoring irradiation again (and tmsting the numerical factors that we just said 
were not to be trusted) we can examine the implications of equation (1221 for pro¬ 
toplanetary disks. Taking Qo = 1.5 and Kq = 2x 10~ 4 cm 2 g- 1 K -2 (5jQj| we find, 
across the region of the disk where ice grains would be the dominant opacity source, 
that, 

r \9/2 


a 


°' 3 (50 Au) 


(123) 


The steep radial dependence of the estimated self-gravitating a means that, unless 
all other sources of transport are extremely small, it will play no role in the in¬ 
ner disk. In the outer disk, on the other hand, we predict vigorous transport. The 
physical origin of the transport is density inhomogeneities that are caused by self¬ 
gravity, which become increasing large as a grows (explicitly, it is found 1861 that 
the average fractional surface density perturbation 8E/E a 1 / 2 ). Since even the 
linear threshold for gravitational instability implies that pressure forces can barely 
resist collapse, we expect that beyond some critical strength of turbulence a self- 
gravitating disk will be unable to maintain a steady-state. Rather, it will fragment 
into bound objects that are not (at least not immediately) subsequently sheared out 
or otherwise disrupted. 

Our discussion up to this point might lead one to conjecture that the threshold for 
fragmentation could be written in terms of a critical dimensionless cooling time. 


fic rit — f- kool ,c 


(124) 


or via a maximum cx cnl that a self-gravitating disk can sustain without fragmenting 
(these are almost equivalent, but defining the threshold in terms of a incorporates 
the varying compressibility as expressed through 7). Gammie (2001) Il2ll . using lo¬ 
cal two dimensional numerical simulations, obtained /3 cl -i t ~ 3 for a two-dimensional 
adiabatic index 7 = 2. Early global simulations B288B . which were broadly con¬ 
sistent with Gammie’s estimate, implied a maximum effective transport efficiency 

acrit-0.1 E9D. 

The idea that the fragmentation threshold is uniquely determined by a single 
number is too simplistic. Several additional physical effects matter. First, fragmenta- 
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tion requires that collapsing clumps can radiate the heat generated by adiabatic com¬ 
pression. There is therefore a dependence not just on the magnitude of the opacity, 
but also on how it scales with density and temperature H66ll85l . Second, if we view 
fragmentation as requiring a critical over-density in a random turbulent field there 
should be a time scale dependence, with statistically rarer fluctuations that lead to 
collapse becoming probable the longer we wait 025 811 . (This introduces an additional 
implicit dependence on y, because the statistics of turbulent density fields depend 
upon how compressible the gas is ED .) Finally, disks can be prompted to fragment 
not only if they cool too quickly, but also if they accrete mass faster than self-gravity 
can transport it away. Accretion-induced fragmentation appears inevitable for very 
massive disks, where it would lead to binary formation 11831 . 

In addition to these physical complexities, there has been debate regarding the 
resolution needed for numerical convergence of simulations of fragmentation 12331 , 
and whether the simplest numerical experiments (which forgo radiation transport 
in lieu of local cooling prescriptions) are even well-posed problems. The numerical 
issues overlap, in part, with the possibility of stochastic fragmentation, and are not 
fully resolved. Disks may fragment for (at least) modestly longer cooling times than 
previously thought, though the high resolution radiation hydrodynamics simulations 
that might definitively resolve the physical questions have yet to be completed. 

Fortunately, if our main interest is in estimating where isolated protoplanetary 
disks ought to fragment, the steep scaling of a(r) (equation 123 i means that we 
can tolerate considerable uncertainty in a cr ; t . For a Solar mass star, fragmentation 
is expected beyond r ~ 10 2 AU f78l 1284 1. with an uncertainty in that estimate of 
perhaps a factor of two. In most (but perhaps not all) cases, it is expected that the 
disk conditions that allow fragmentation would lead to objects with masses in the 
brown dwarf regime, or above 1 1841 . 


4.3 Magnetohydrodynamic turbulence and transport 

The Rayleigh stability criterion (equation |110| > applies to a fluid disk. It does not 
apply to a disk containing even an arbitrarily weak magnetic field, if that field is 
perfectly coupled to the gas (the regime of ideal MHD). In ideal MHD, indeed, a 
weakly magnetized disk has entirely different stability properties from an unmag¬ 
netized one, and is unstable provided only that the angular velocity decreases out¬ 
ward. This is the magnetorotational instability (MRI) |37][38l, which is accepted as 
the dominant source of turbulence in well-ionized accretion disks (winds could still 
contribute to or dominate angular momentum loss). In protoplanetary disks the ideal 
MHD version of the MRI applies only in the thermally ionized region close to the 
star; across most of the disk we also need to consider both the dissipative (Ohmic 
diffusion, ambipolar diffusion) and the non-dissipative (the Hall effect) effects of 
non-ideal MHD. The Ohmic and ambipolar terms can be considered as modifying 
— albeit very dramatically — the ideal MHD MRI, while the Hall term introduces 
new effects (in part) via the Hall shear instability 11871 , which is a different beast 
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Fig. 17 Illustration showing why a weak vertical magnetic field destabilizes a Keplerian disk (the 
magnetorotational instability 1 1381 ). An initially uniform vertical field (weak enough that magnetic 
tension is not dominant) is perturbed radially. Due to the shear in the disk, an inner fluid element 
coupled to the field advances azimuthally faster than an outer one. Magnetic tension along the field 
line then acts to remove angular momentum from the inner element, and add angular momentum 
to the outer one. This causes further radial displacement, leading to an instability. 


unrelated to the ideal MHD MRI. The phenomenology of disk instabilities in non¬ 
ideal MHD is rich, and appears to give rise to both turbulent and laminar angular 
momentum transport as well as phenomena, such as MHD disk winds, that may be 
observable. 


4.4 The magnetorotational instability 


The MRI El is an instability of cylindrical shear flows that contain a weak 
(roughly, if the field is vertical, sub-thermal) magnetic fielcj^] In ideal MHD the 
condition for instability is simply that, 


df2 2 

~dr~ 


< 0 . 


( 125 ) 


The fact that this condition is always satisfied in disks (though not in star-disk 
boundary layers) accounts for the MRI’s central role in modern accretion theory. 


6 The mathematics of the MRI was worked out by Velikhov 134H and Chandrasekhar EO around 
1960. Thirty years passed before Balbus & Hawley 03 recognized the importance of the instabil¬ 
ity for accretion flows. 
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Figure [TT] illustrates what is going on to destabilize a disk that contains a mag¬ 
netic field. A basically vertical field is slightly perturbed radially, so that it links 
fluid elements in the disk at different radii. Because of the shear in the disk, the 
fluid closer to the star orbits faster than the fluid further out, creating a toroidal field 
component out of what was initially just vertical and radial field. The tension in the 
magnetic field linking the two elements (which can be thought of, even mathemati¬ 
cally, as being analogous to a stretched spring) imparts azimuthal forces to both the 
inner fluid (in the direction opposite to its orbital motion) and the outer fluid (along 
its orbital motion). The tension force thus reduces the angular momentum of the in¬ 
ner fluid element, and increases that of the outer element. The inner fluid then moves 
further inward (and the outer fluid further outward) and we have an instability. 

We can derive the MRI instability condition in a very similar setup as Figure o 
Consider a disk with a power-law angular velocity profile, Q « r q , that is threaded 
by a uniform vertical magnetic field Bq. We ignore any radial or vertical variation 
in density (and consistent with that, ignore the vertical component of gravity) and 
adopt an isothermal equation of state, P = pc 2 , with c s a constant. Our task is to 
determine whether infinitesimal perturbations to this equilibrium state are stable, or 
whether instead they grow exponentially with time, signaling a linear instability. 

To proceed (largely following 033 ) we define a locally Cartesian patch of disk 
that corotates at radius ro, where the angular frequency is Qq. The Cartesian co¬ 
ordinates (x.y.z) are related to cylindrical co-ordinates (r, @.z r ) via, 


* = r— ro, 

y = r 0 </>, 

z = z- (126) 


The local “shearing-sheet” (or in three dimensions, “shearing box”) model is useful 
for both analytic stability studies, and for numerical simulations. In this co-rotating 
frame, the equations of ideal MHD pick up terms representing the fictitious Coriolis 
and centrifugal forces, 


f + V.(pv) 
A +(V .V), 

dB 

dt 


0 , 

1 1 , 

- VP -\ -(V x B) x B — 2Qo x v + IqQfixx, 

p 4 Tip 

Vx(vxB). (127) 


Here x is a unit vector in the x-direction. As noted above, the initial equilibrium has 
uniform density, p = po, and a magnetic field B = (0,0,Bo)- There are no pressure 
or magnetic forces, so the velocity field is determined by a balance between the 
Coriolis and centrifugal terms. 


2Qq x v = 2qQ. 2 ) xk. 


(128) 
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The equilibrium velocity field that completes the definition of the initial state is, 

v= (0,-gApqO), (129) 

which has a linear shear (with q = 3/2 for a Keplerian disk) around the reference 
radius ro- 

To assess the stability of the equilibrium, we write the density, velocity and mag¬ 
netic field as the sum of their equilibrium values plus a perturbation. We can recover 
the MRI with a particularly simple perturbation which depends on z and t onlyj^] For 
the velocity components, for example, we write, 

V.v = v' x {z,t), 

Vy = -qQox + v' y (z,t), 

v z = v’ z (z,t), (130) 

and do likewise for the density and magnetic field. We substitute these expressions 
into the continuity, momentum and induction equations, and discard any terms that 
are quadratic in the primed variables, assuming them to be small perturbations. This 
would give us seven equations in total (one from the continuity equation, and three 
each from the other equations), but the x and y components of the momentum and 
induction equations are all we need to derive the MRI. The relevant linearized equa¬ 
tions are. 


dt 


dv' x 

B 0 dB' x 

dt 

47rp 0 dz 


Bo dB' y 

q£2ov' x 

4Kpo dz 

dB' 

dv' x 

dt 


dB[. 

dv[, 

dt 

= Bo J7-1 


+ 2Qq vl 


3 ” 


— 2£2qv', 


(131) 


We convert these linearized differential equations into algebraic equations by taking 
the perturbations to have the form, e.g., 

B' x = B' x e^ m ~ kz \ (132) 


where ft) is the frequency of a perturbation with vertical wave-number k. The time 
derivatives then pull down a factor of /ft), while the spatial derivatives become ik. 
Our four equations simplify to, 


7 An analysis that retains the ^-dependence can be found in the original Balbus & Hawley 
(1991) paper (37], and follows an essentially identical approach. Studying the stability of non- 
axisymmetric perturbations (in y), however, requires a different and more involved analysis 
(87 ll25Tll326l[267l . 
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. BqB' . 

i(Ov' = — ik -- + 2Qqv, 

47rp 0 v 

, B 0 B' 

iav = -ik—— + (q- 2)£2 0 v x , 

47ip 0 

icoB' x = - ikBo v r x , 

i( 0 B' y =-ikB 0 V y -qn 0 B' x . (133) 


(We’ve dropped the bars on the variables for clarity.) Eliminating the perturbation 
variables from these equations, we finally obtain the MRI dispersion relation, 

ft) 4 — (0 2 [2k 2 v 2 A + 2(2-q)D.l] +k 2 v\ [k 2 v 2 A -2qQl] = 0, (134) 


where v A = B 2 J (47rpo) is the Alfven speed associated with the net field. 

If ft) 2 > 0 then co itself will be real and the perturbation e ,cot will oscillate in 
time. Instability requires ft) 2 < 0, since in this case ft) is imaginary and the perturba¬ 
tion will grow exponentially. Solving the dispersion relation we find the instability 
criterion is, 


(kv A ) — 2pf2 ( f < 0. 


(135) 


Letting the field strength go to zero (B- —> 0, v A —> 0) we find that the condition for 
instability is simply that q > 0, i.e. that the angular velocity decrease outward. Even 
for an arbitrarily weak field, the result is completely different from Rayleigh’s for a 
strictly hydrodynamic disk. 

The growth rate of the instability and what it means for the magnetic field to 
be “weak” can also be derived from equation ( 1 34[ ). Specializing to a Keplerian 
rotation law with q = 3/2 the dispersion relation takes the form shown in Figure [j~8| 
For a fixed magnetic field strength (and hence a fixed Alfven speed v A ) the flow is 
unstable for wavenumbers k < k a ll (i.e. on large enough spatial scales), where, 


^crit l ’A — 5/3 f2o- 


(136) 


As the magnetic field becomes stronger, the smallest scale A = 2n/k cAx which is un¬ 
stable grows, until eventually it exceeds the disk’s vertical extent « 2h. For stronger 
vertical fields no unstable MRI modes fit within the disk, and the instability is sup¬ 
pressed. Using h = c s /Q, the condition that the vertical magnetic field is weak 
enough to admit the MRI (i.e. that A < 2 h) becomes 

Bl < —pc 2 (137) 

K 

If we define the plasma j3 parameter as the ratio of gas to magnetic pressure. 




8 7tP 
B 2 


( 138 ) 


this condition can be expressed alternatively as, 
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Fig. 18 The unstable branch of the MRI dispersion relation is plotted for a Keplerian rotation law. 
The flow is unstable (ar < 0) for all spatial scales smaller than Icva < y/3£2 (rightmost dashed 
vertical line). The most unstable scale (shown as the dashed vertical line at the center of the plot) 
is close to Icva — Q ■ 


J3> 


2 n 2 

T~ 


(139) 


A magnetic field whose vertical component approaches equipartition with the ther¬ 
mal pressure (/3 ~ 1) will be too strong to admit the existence of linear MRI modes, 
but a wide range of weaker fields are acceptable. 

The maximum growth rate is determined by setting dot) 2 /d(kv a) = 0 for the un¬ 
stable branch of the dispersion relation plotted in figure [18] The most unstable scale 
for a Keplerian disk is. 


(Mm ax 



(140) 


where the growth rate is, 


| tAnax | — ^ ^0 • 


(141) 


This result implies an extremely vigorous growth of the instability, with an exponen¬ 
tial growth time scale that is a fraction of an orbital period. This means that if a disk 
is unstable to the MRI it will invariably dominate the evolution. 
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4.4.1 Non-ideal MHD 


The MRI in its ideal MHD guise is rele vant to protoplanetary disks only in the ther¬ 
mally ionized region close to the star (f 2.3.11, where T > 10 3 K. The very weakly 
ionized gas further out is imperfectly coupled to the magnetic field, and this both 
modifies the properties of the MRI and leads to new MHD instabilities. We will 
begin by sketching the derivation of the non-ideal MHD equations (following Bal- 
bus li36l . who justifies several of the approximations that we will make), and then 
estimate the magnitude of the extra terms that arise in protoplanetary disks. 

The physics of how magnetic fields affect weakly-ionized fluids is easy to visu¬ 
alize. We consider a gas that is almost entirely neutral, with only a small admixture 
of ions and electrons (analogous considerations apply if the charge carriers are dust 
particles, but we will not go there). Magnetic fields exert Lorentz forces on the 
charged species, but not on the neutrals. Collisions between the neutrals and either 
the ions or the electrons lead to momentum exchange whenever the neutral fluid has 
a velocity differential with respect to the charged fluids. 

We begin by considering the momentum equation. For the neutrals we have. 


d\ 

p-^+p(v-v)v = -'yp-pv®- PnI -p ne . 


(142) 


Here p, v and P (without subscripts) refer to the neutral fluid, and /;„/ and p ne are 
the rate of momentum exchange due to collisions between the neutrals and the ions 
/ electrons respectively. Identical equations apply to the charged species, but for the 
addition of Lorentz forces. 


P,- S f+P,(y,-V)y, 

dvr , 

P/-^-+P/(v/-V)v/ 


—VRe — p e V <J> — en e + Vg * B ^ - p e „, 

— VP/ — p/V<J> +Zen/ ^E+ —j —pin- (143) 


In these equations E and B are the electric and magnetic fields, the ions have charge 
Ze, where — e is the charge on an electron, and of course p ne = —p en and p„/ = —pi„. 
Having three momentum equations looks complicated, but we can make a large 
simplification to the system by noting that the time scale for macroscopic evolution 
of the fluid is generally much longer than the time scale for collisional or magnetic 
forces to alter a charged particle’s momentum. We can then ignore everything in the 
charged species’ momentum equations, except for the Lorentz and collisional terms. 
For the ions we have, 


Zen, (^E+ V ^ B j - Pr „ = Q, (144) 

with a similar equation for the electrons. Imposing charge neutrality, n e = Zni, we 
eliminate the electric field between the ion and electron equations to find an expres¬ 
sion for the sum of the momentum transfer terms. 
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Pin + Pen = - (V/ - \ e ) X B. 

C 

The current density J = en e (\j — \ e ), so we can write this as, 

J xB 

Pin T Pen — 


(145) 


(146) 


Finally, we go to Maxwell’s equations, and note that the current can be written as, 

(147) 


An „ „ 1 <9E 

—J = V x BH-=—. 

c c at 


The second term in Maxwell’s equation is the displacement current, which is 
ff{v 2 /c 1 ) and consistently ignorable in non-relativistic MHD. Doing so, we sub¬ 
stitute equation (|146|l in the neutral equation of motion to obtain. 


<9v 1 

p — +p(v-V)v = -VP-pV<£ + — (V x B) x B. 
at An 


(148) 


This is identical to the ideal MHD momentum equation (stated without derivation as 
equation 1 127| > and pleasingly simple; we have reduced the three momentum equa¬ 
tions to an equation for a single (neutral) fluid with a magnetic force term whose 
dependence on B is independent of the make-up of the gas. All of the complexities 
enter only via the induction equation. 


The consistent simplification of non-ideal MHD to a momentum equation for 
a single fluid is not always possible. Roughly speaking it works provided that 
the plasma’s inertia is negligible compared to that of the neutral fluid, the 
coupling between charged and neutral species is strong, and the recombina¬ 
tion time is short. Zweibel [1364ij gives an accessible account of the conditions 
necessary for a valid single-fluid description. Although some of the early an¬ 
alytic and numerical work on the MRI in weakly-ionized disks utilized a two- 
fluid approach f56lll4H , in many protoplanetary disk situations a single fluid 
model is both justified 1 1251 and substantially simpler. 


Deducing the non-ideal induction equation requires us to specify the form of the 
momentum coupling terms. Writing these in standard notation (which is different 
for the two terms, somewhat obscuring the symmetry), 

Pne = n e v ne m e (\ V,.), 

Pm = PP/7(v-V/), (149) 

where v ne is the collision frequency of an electron with the neutrals, and y is called 
the drag-coefficient. The ion-neutral coupling involves longer-range interactions 
than the electron-neutral coupling, and is accordingly stronger ll36ll . 
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We now go back to the force balance deduced from the electron momentum equa¬ 
tion, 

- en e + Vg * B ^ - Pen = 0, (150) 

and attempt to write the terms involving y e and p en entirely in terms of the current. 
We start with the exactly equivalent expression, 

I v/ vn 

E +-[v + (v e -v/) + (v/-v)] xB+-y^[(y e -v/) + (v/-v)] =0, (151) 

and deal with the terms in turn. We have two terms that involve (v e — V/), which can 
be replaced immediately with the current. 


(v e -v 7 ) 


,1 

en e 


(152) 


The first term with (v/ — v) is exactly equal to pj n / ( ppij )■ If, however, \pj n l»l P en 1 9 
then equation (| 146|> implies that, approximately. 


(v/-v) 


JxB 

cppij' 


(153) 


Finally it can be shown (see Balbus f36l for details) that the final term with (v/ — v) 
can be consistently dropped. The version of Ohm’s Law that we end up with is, 


vxB JxB (JxB)xB 

E H-f ■ 




en e c 


c 2 PPiY 


-3 = 0. 


(154) 


The non-ideal induction equation is then obtained by applying Faraday’s law, 


V x E = — 


1 <9B 
c dt ’ 


(155) 


to eliminate any explicit reference to the electric field. In its usual form. 


c)\\ „ 

= Vx 

dt 


v x B r/V x B 


JxB (J x B) x B 

en e cypp! 


(156) 


We have defined the magnetic resistivity, 

c 2 

0 =- 

4kg 

where a is (here) the electrical conductivity, 

e 2 n e 

G = -. 

Hie ^en 

As before we can replace the current with the magnetic field via. 


057) 


058) 
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J=fVxB, (159) 

4 n 

so that the induction equation is solely a function of B. The terms on the right-hand- 
side are referred to as the inductive. Ohmic, Hall and ambipolar terms respectively. 

The non-ideal terms in the induction equation depend upon the ionization state 
of the gas (through n e and p/) and upon the collision rates between the neutral and 
charged species (via T) and y). Standard values for these quantities are l56l i95;l. 

T 7 = 234 ^ — ^ T 1/2 cm 2 s _1 , 

y = 3x 10 13 cm 3 s _I g _1 . (160) 

We are now ready to estimate the importance of the non-ideal terms in the environ¬ 
ment of protoplanetary disks, and to ask what effect they have both on the MRI, and 
on the more general question of whether there is MHD turbulence or transport in 
disks. 


4.4.2 Ohmic, ambipolar and Hall physics in protoplanetary disks 


The non-ideal terms in equation (156 1 all depend inversely on the electron or ion 
density, so the strength of all non-ideal MHD effects relative to the inductive term 
increases with smaller ionization fraction. The three terms also have different depen¬ 
dencies on density, magnetic field strength and temperature, so the relative ordering 
of the non-ideal MHD effects varies with these parameters. 

The Ohmic, Hall and ambipolar terms have different dependencies on the mag¬ 
netic field geometry, and in a disk setting they influence the MRI in distinct ways 
(most importantly, the Hall effect differs from the others in being non-dissipative). 
There is therefore no model-independent way to precisely demarcate when each 
term will affect disk evolution. As a first guess, however, we can treat the magnetic 
field as a scalar and simply take the ratio of the Hall to the Ohmic term and the 
ambipolar to the Hall term. 


H cB 
O Alter) n e ' 
A en e B 

H cypp/' 


061) 


Since rj (n/n e ) and n e p/ both of these ratios depend on ( B/n ). Substituting for 
V) and y, and taking the ion mass that enters into the ambipolar term as 30 inn, we 
can estimate the magnetic field strength for which the Ohmic and Hall terms have 
equal magnitude, and similarly for the Hall and ambipolar terms. 
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n / cm -3 


Fig. 19 Regions of the (n,B) parameter space in which different non-ideal terms are dominant. 
The boundary between the Ohmic and Hall regimes is plotted for T = 100 K (solid line) and also 
for temperatures of 10 3 K (upper dashed line) and 10 K (lower dashed line). The red line shows 
a very rough estimate of how the magnetic field in the disk might vary with density between the 
inner disk at 0.1 AU and the outer disk at 100 AU. 

b »-** 4x10 " 3 (toi 4^) g ' (i62) 

Figure[l9]shows these dividing lines in the (n. B) plane. Ohmic diffusion is dominant 
at high densities / low magnetic field strengths. Ambipolar diffusion dominates for 
low densities / high field strengths. The Hall effect is strongest for a fairly broad 
range of intermediate densities. 

Estimating where protoplanetary disks fall in the (n,B) plane can be done in 
various ways. For a Solar System-motivated estimate we can start with the disk field 
inferred from laboratory measurements of chondrules in the Semarkona meteorite 
cm, which suggest that near the snow line (r ~ 3 AU) the disk field was B ~ 
0.5 G. (There are caveats and a large systematic uncertainty associated with this 
measurement, all of which we ignore for now.) Let us assume that the surface density 
and temperature profiles are Z ^ r -3/ 2 and T r 1 /2 respectively, and that the 
magnetic field pressure is the same fraction of the thermal pressure at all radii in 
the disk. Taking E ~ 300 g cm ~ 2 and T = 150 K at 3 AU, the inferred scalings of 
mid-plane density and magnetic field strength are then. 



( 163 ) 
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The track defined by these relations is plotted in Figure[l9]for radii between 0.1 AU 
and 100 AU. 

As we have emphasized, neither our approach to ranking the strength of the non¬ 
ideal terms, nor our estimate of the radial scaling of disk conditions, are anything 
more than crude guesses. Other approximations are equally valid (for example, one 
can order the terms in the (n,T) plane instead ll4T[|l88[|T6l ). Nevertheless, because 
n and B vary by so many orders of magnitude across Figure[l9]the critical inferences 
we can draw are quite robust. We predict that the Hall effect is the dominant non¬ 
ideal MHD process at the disk mid-plane between (conservatively) 1 AU and 10 AU. 
Ohmic diffusion can become important as we approach the thermally ionized region 
interior to 1 AU. Ambipolar diffusion dominates at sufficiently large radii, of the 
order of 100 AU, and in the lower density gas away from the mid-plane. 


4.4.3 The dead zone 


The linear stability of Keplerian disk flow in non-ideal MHD has been extensively 
investigated (see, e.g. Il56l 1159113441 l4Tl 188 ( JT|). and the reader interested in the 


non-ideal analogs to the MRI dispersion relation derived as equation (134 1 should 


start there. Proceeding less formally, we follow Gammie (1996) HI 191 to estimate the 
conditions under which Ohmic dissipation (ignoring for now the Hall term) would 
damp the MRI. The basic idea is to compare the time scale on which the ideal MRI 
would generate tangled magnetic fields to that on which Ohmic diffusion would 
smooth them out. We first note that diffusion erases small-scale structure in the 
field more efficiently than large scale features, so that the appropriate comparison 
is between growth and damping of the largest scale MRI models. Starting from the 
MRI dispersion relation (equation |134| i, for a Keplerian disk we consider the weak- 
field / long wavelength limit (kv,\ j Cl <C 1). The growth rate of the MRI is. 


|co| ~ x/3 kv A . 


(164) 


Writing this as a function of the spatial scale A = 2n/k, we have, 

\co\~2nV3j. (165) 

Up to numerical factors the MRI on a given scale then grows on the Alfven crossing 
time. Equating this growth rate to the Ohmic damping rate. 


I ®j? I ~ 


A 2 


(166) 


we conclude that Ohmic dissipation will suppress the MRI on the largest available 
scale A « h provided that, 


T) > 2nV2>VAh. 


(167) 
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We can express this result in a different form. Analogous to the fluid Reynolds 
number (equation [108]) the magnetic Reynolds number Re v/ is defined as, 


Re*/ = 


UL 

V 


(168) 


where U is a characteristic velocity and L a characteristic scale. Taking U = va and 
L = h for a disk, the condition for Ohmic dissipation to suppress the MRI becomes 


Re*/ < 1 


069) 


where order unity numerical factors have been omitted. 

We now convert the condition for the suppression of the MRI into a limit on the 
ionization fraction x = n e ln. We make use of the formula for the magnetic resistiv¬ 
ity (equation 160 1 and assume that Maxwell stresses transport angu lar momentum, 
so that a ~ /cj (this follows approximately from equation 109 1 . The magnetic 
Reynolds number can then be estimated to be. 


Re M — 


VAh 

V 


a 1/,2 c 2 

rjQ 


(170) 


Substituting for 7] and c 2 , the magnetic Reynolds number in a protoplanetary disk 
scales as. 


Rev/ 


1.4 x 10 12 x 



1/2 



T \ 1/2 /M*\ 1/2 
300 Kj \M^J 


(171) 


For the given parameters, the critical ionization fraction below which Ohmic diffu¬ 
sion will quench the MRI is 

Xcrit-10- 12 . (172) 

Clearly a very small ionization fraction suffices to couple the magnetic field to the 
gas and allows the MRI to operate, but there are large regions of the disk where even 
these ionization levels are not obtained and non-ideal effects are important. 

Based on this analysis Gammie 111 191 noted, first, that the criterion for the MRI 
to operate under near-ideal MHD conditions in the inner disk coincides with the re¬ 
quirement that the alkali metals are thermally ionized (Figure [6] . The development 
of magnetized turbulence in the disk at radii where T > 10 3 K can therefore be mod¬ 
eled in ideal MHD. Second, he proposed that Ohmic diffusion would damp MHD 
turbulence in the low ionization environment near the disk mid-plane on scales of 
the order of 1 AU, creating a dead zone of sharply reduced turbulence and transport. 
Gammie’s original model is incomplete, as it did not include either ambipolar dif¬ 
fusion or the Hall effect, but the basic idea motivates much of the current work on 
MHD instabilities in disks. 
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4.4.4 Turbulence and transport in non-ideal MHD 


Once we consider the full set on non-ideal terms, the first question is to assess the 
level of turbulence and transport that is expected as a function of their strengths. 
This is a well-defined but already difficult theoretical question, given that interesting 
values of the Ohmic, ambipolar and Hall terms span a broad range (depending phys¬ 
ically on the temperature, density, ionization fraction and magnetic field strength). 
To define the problem in its most idealized form, we rewrite the non-ideal induction 
equation ( 156} as. 


cJK ^ 

-2~= VX 

dt 


JxB (JxB)xB 

vx B — r]o V xB-i)// —-- rj A 


B 


B 2 


(173) 


where dimensionally T)q, 7]// and r\ A are all diffusivities. There are different ways to 
construct dimensionless numbers from the diffusivities, but one useful set is. 


Aq = 


Qrio 


(Ohmic Elsasser number), 


Ha = 


V A 


Qt]h 


(Hall Elsasser number), 


Am = 


'A 


QriA 


(ambipolar Elsasser number). 


(174) 


(Note that the ambipolar Elsasser number can also be written as Am = jpi/Q, 
which is the number of ion-neutral collisions per dynamical time .) We further 
specify the net vertical magnetic field (if any) via the ratio of the mid-plane gas and 
magnetic field pressures (equation|138|l. 


Pz = 


%nP 

~w 


(175) 


Our question can then be rephrased; what is the level of angular momentum trans¬ 
port and turbulence in an MHD disk at radii where the non-ideal terms are charac¬ 
terized by the dimensionless parameters (Aq, Ha, Am, fi z ). 

Ohmic diffusion acts as a strictly dissipative process that stabilizes disks to mag¬ 
netic field instabilities on scales below some critical value. Ambipolar diffusion is 
in principle more complex, because it does not dissipate currents that are paral¬ 
lel to the magnetic field. This distinction substantially impacts the linear stability 
of ambipolar-dominated disks II1881 . but appears to matter less for the non-linear 
evolution, whose properties are analogous to Ohmic diffusion. For both dissipative 
processes, simulations show that MRI-driven turbulence is strongly damped when 
the relevant dimensionless parameter (either Ao or Am) drops below a critical value 
that depends upon the initial field geometry but is ~ 1 — 10 2 l30QII335ll310ll32ll308L 
Consistent with the dead zone idea 111 1911 . we therefore expect substantial modifica¬ 
tion of the properties of MHD turbulence both in the mid-plane around 1 AU (where 
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Fig. 20 An illustration (adapted from Lesur, unpublished) of the Hall shear instability 11871 . In the 
presence of a vertical field threading the disk, the Hall effect acts to rotate an initially toroidal field 
component either clockwise or anti-clockwise, depending upon the sign of the vertical field. The 
rotated field vector is then either amplified or damped by the shear. This instability differs from 
the MRI both technically (in that there is no reference to orbital motion, any shear flow suffices) 
and physically (via the dependence on the direction of the vertical magnetic field as well as its 
strength). 


Ohmic diffusion is the dominant dissipative process) and in the disk atmosphere and 
at large radii ~ 10 2 AU (where ambipolar diffusion dominates). 

The MRI dispersion relation is also modified by the Hall effect [1 34411411 . which 
differs from the other non-ideal terms in that it modifies the field structure without 
any attendant dissipation. In this respect the Hall term most closely resembles the 
inductive term V x (v x B), and its strength can usefully be characterized by the 
ratio of the Hall to the inductive term. Non-linear simulations of the Hall effect in 
disks, which were pioneered by Sano & Stone 13011, 1302 1. have only recently been 
able to access the strongly Hall-dominated regime relevant to protoplanetary disks 
1118911195112711281 . In vertically stratified disks with a net vertical field, Lesur et al. 
(2014) 11951 find that the Hall effect has a controlling influence on disk dynamics 
on scales between 1 AU and 10 AU. For /3- = 10 5 a strong but laminar Maxwell 
stress (i.e. one dominated by large-scale radial and toroidal fields in equation 109 1 
results when the net field is aligned with the rotation axis of the disk, whereas anti¬ 
alignment leads to extremely weak turbulence and transport. 

The results from Hall-MHD simulations of protoplanetary disks are best inter¬ 
preted not as a modification of the MRI, but rather as the signature of a distinct 
Hall shear instability mm In the presence of a net vertical magnetic field the Hall 
effect acts to rotate magnetic field vectors lying in the orbital plane (Figure [20]), 
with the sense of the rotation determined by whether the new held is aligned or 
anti-aligned to the rotation axis. In the aligned held case, the Hall-induced rotation 
allows the magnetic held to be amplihed by the shear, while damping occurs in the 
anti-aligned limit. Unlike the MRI, the Hall shear instability does not depend on 
the Coriolis force, and is indifferent to the sign of the angular velocity gradient. 
By generating a radial held directly from an azimuthal one, the Hall effect (given a 






























66 


Philip J. Armitage 


net field) supports a mean-field disk dynamo cycle 1 329111951 that is qualitatively 
different from anything that is possible in ideal MHD. 

Elementary arguments show that the Hall effect is important across a range of 
radii in protoplanetary disks (Figure [T9}. Whether existing simulations cap¬ 
ture the full extent of non-ideal MHD behavior in disks is, however, open to 
question. The numerical implementation of the Hall effect in simulation codes 
poses significant challenges, and the presence of large-scale fields in the satu¬ 
rated state suggests that local simulations may not be adequate to describe the 
outcome. The level of turbulence that accompanies the predominantly large- 
scale transport by Maxwell stresses appears particularly uncertain. 


Going beyond the idealized question of the effect of the non-ideal terms on tur¬ 
bulence and transport, our goal is to use the results described above to predict the 
structure and evolution of protoplanetary disks. This introduces new layers of un¬ 
certainty. To predict disk properties from first principles, we need at a minimum 
to know the strength of the different sources of ionization ( §2.3.2 1 , the rates of gas- 
phase and dust-induced recombination, and the global evolution of any net magnetic 
held ( §3.5. 1 [ ). We also need to model (or have good reasons to ignore) various non- 
MHD effects, including the hydrodynamic angular momentum transport processes 
already discussed and mass loss by photoevaporation [0 (Sj8}. Given these uncer¬ 
tainties, the best that we can currently do is to highlight a number of qualitative 
predictions that receive support from numerical simulations: 


• Net vertical magnetic fields are important for disk evolution. A vertical mag¬ 
netic held enhances angular momentum transport by the MRI even in ideal MHD 

1 jr\ 

(roughly as a °c fj 7 ~ 111401 ). In non-ideal MHD, local simulations by Simon et 
ah 1 30811307 ] suggest that ambipolar damping of the MRI in the outer disk pre¬ 
vents resupply of the inner disk with gas unless a net held is presenj^] A net held 
with J 3 Z ~ 10 4 — 10 5 suffices, comparable to the helds expected in global models 
of flux evolution 1134ft but much weaker than the likely initial held left over from 
star formation. Accretion on scales of 30-100 AU occurs predominantly through 
a thin surface layer that is ionized by FUV photons 1266,308 3071. and is largely 
independent of the Hall effect l28l . 

• MHD winds and viscous transport can co-exist in disks. Local numerical simu¬ 
lations in ideal MHD by Suzuki & Inutsuka 13201 showed that in the presence of 
net vertical held, the MRI was accompanied by mass and angular momentum loss 
in a disk wind. Winds are likewise seen in net held protoplanetary disk simula¬ 
tions at 1 AU that include Ohmic and ambipolar diffusion ll33l . in simulations at 


8 It is not obvious that the inner disk is resupplied by gas, or, to put it more formally, that the disk 
attains a steady-state. Out to ~ 10 AU the viscous time scale is short enough that the disk will 
plausibly adjust to a steady state (provided only that a steady state is possible, see §5j, but no such 
argument works out to 100 AU. Ultimately the question of whether gas at 100 AU ever reaches the 
star will need to be settled by observations as well as by theory. 
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30-100 AU where ambipolar diffusion dominates 13071 . and in simulations that 
include all non-ideal terms urn Caution is required before interpreting these 
local simulation results as quantitative predictions, because although the effec¬ 
tive potential for wind launching is correctly represented (equation |97j) there is 
a known and unphysical dependence of the mass loss rate on the vertical size of 
the simulation domain C33. Outflows are also seen in the (few) existing global 
net field simulations 11311 . however, supporting the view that weak field MHD 
outflows could be a generic feature of protoplanetary disk accretion. Drawing on 
these results, Bai l26l has proposed that MHD winds may dominate disk evolu¬ 
tion on AU-scales, and co-exist with a mixture of viscous transport and wind loss 
on larger scales in the ambipolar dominated regime 113071 . 

• Turbulence and angular momentum transport are not synonymous. In classical 
disk theory, the value of a determines not only the rate at which the disk evolve, 
but also the strength of turbulence and its effect on small solid particles. This 
link is doubly broken in more complete disk models. First, as already noted, an¬ 
gular momentum loss via winds (which need not be accompanied by turbulence) 
may be stronger than viscous transport at some radii. Second, even the internal 
component of transport may be primarily a large-scale “laminar” Maxwell stress, 
rather than small-scale turbulence 11307111951 . 

• The sign of the net field could lead to bimodality in disk properties. The Hall 
effect is the strongest non-ideal term interior to about 10 AU, and simulations 
111951 l27l 28) confirm the expectation from linear theory 1 3441 [4X111871 that a 
disk with a weak field that is aligned to the rotation axis behaves quite differently 
from one with an anti-aligned field. Although there are possible confounding 
factors — for example the long-term evolution of the net field may itself differ 
with the sign of the field — it appears likely that the striking asymmetry seen in 
simulations introduces some observable bimodality in disk structure 13111 . 

Figure [21] illustrates a possible disk structure implied by the above results. The 
figure should be regarded as a work in progress; there is plenty of work remain¬ 
ing before we fully understand either the physics of potential angular momentum 
transport and loss processes, or how to tie that knowledge together into a consistent 
scenario for disk structure and evolution. 


Hall MHD can also affect the collapse of molecular clouds and the formation 
of protostellar disks 1182113301 . influencing for example their initial sizes. 
The Hall current, along with the other non-ideal terms, can also modify the 
accretion properties of circumplanetary disks 111701 
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0.1 AU 1 AU 10AU 100 AU 

Fig. 21 A suggested structure for protoplanetary disks if MHD processes dominate over other 
sources of transport. The different regions are defined by the strength of non-ideal MHD terms 
(Ohmic diffusion, ambipolar diffusion and the Hall effect), and by mass and angular momentum 
loss in MHD disk winds. The Hall effect is predicted to behave differently if the net field threading 
the disk is anti-aligned to the rotation axis (here, alignment is assumed). Ionization by stellar X- 
rays and by FUV photons couples the stellar properties to those of the disk. 


4.5 Transport in the boundary layer 

The nature of transport in the boundary layer deserves a brief discussion. As dis¬ 
cussed in j p.2.1| boundary layers are expected when the accretion rate is high enough 
to overwhelm the disruptive influence of the stellar magnetosphere (see equation [75] 
for a semi-quantitative statement of this condition). For most stars this requires high 
accretion rates, so the boundary layer and adjacent disk are hot enough to put us 
into the regime of thermal ionization and ideal MHD. In the disk, we then expect 
angular momentum transport via the MRI. In the boundary layer, however, we face 
a problem. By definition, dfl/'dr > 0 in the boundary layer (see Figure [9ji, and this 
angular velocity profile is stable against the MRI (equation |1 25 | l. Something else 
much be responsible for transport in this region. 

The angular velocity profile in the boundary layer is stable against the genera¬ 
tion of turbulence by either the Rayleigh criterion or the MRI. It turns out to be 
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unstable, however, to the generation of waves via a mechanism analogous to the 
hydrodynamic Papaloizou-Pringle instability of narrow tori 12641 . Belyaev et al. 
Eum, using both analytic and numerical arguments, have shown that waves gen¬ 
erated from the supersonic shear provide non-local transport of angular momentum 
(and energy) across the boundary layer. Magnetic fields are amplified by the shear 
IT3l but do not play an essential role in boundary layer transport l52l . In protostellar 
systems boundary layers are present during eruptive accretion phases (see |5]l when 
strong radiation fields are present fl76l . Future work will need to combine the recent 
appreciation of the importance of wave angular momentum transport with radiation 
hydrodynamics for a full description of the boundary layer. 


5 Episodic accretion 

Young stellar objects (YSOs) are observed to be variable. The short time scale (last¬ 
ing hours to weeks) component of that variability is complex |84| . but can probably 
be attributed to a combination of turbulent inhomogeneities in the inner disk, stellar 
rotation ED- and the complex dynamics of magnetospheric accretion El- There is 
also longer time scale variability — lasting from years to (at least) many decades, 
that in some cases takes the form of well-defined outbursts in which the YSO bright¬ 
ens dramatically. The traditional classification of outbursting sources divides them 
into FU Orionis events ms 1 139 1. characterized by a brightening of typically 5 
magnitudes followed by a decay over decades, and EXors US), which display re¬ 
peated brightenings of several magnitudes over shorter time scales. The statistics on 
these uncommon long-duration outbursts (especially FU Orionis events 1114811 ) are 
limited, and it is not even clear — either observationally or theoretically — whether 
FUOrs and EXors are variations on a theme or genuinely different phenomena ETI . 
Nonetheless, it is established that episodic accretion is common enough to matter 
for both stellar accretion and for planet formation processes occurring in the inner 
disk ED. Our focus here is on the origin of these accretion outbursts. 

Observations show that FU Orionis outbursts involve a large increase in the mass 
accretion rate through the inner disk on to the star 11391 . During the outburst the 
inner disk will be relatively thick (h/r ~ 0.1), and hot enough to be thermally ion¬ 
ized. We therefore expect efficient angular momentum transport from the MRI, with 
a ss 0.02. Writing the viscous time scale (equation[58j) in terms of these parameters. 



(176) 


we can estimate the disk radius associated with a (viscously driven) outburst of 
duration tburst. 



(177) 


For a Solar mass star we find, 
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Fig. 22 An illustration of some of the processes suggested as the origin of episodic YSO accretion. 



Disk-driven outbursts of broadly the right duration are thus likely to involve events 
on sub-AU scales, and could be associated physically with the magnetosphere, with 
the thermally ionized inner disk, or with the inner edge of the dead zone. 

The physical origin of episodic accretion in YSOs has not been securely identi¬ 
fied. Mooted ideas, illustrated in Figure [22] fall into two categories. The first cate¬ 
gory invokes secular instabilities of protoplanetary disk structure that may occur on 
AU and sub-AU scales. The idea is that the inner disk may be intrinsically unable 
to accrete at a steady rate, and instead alternates between periods of high accretion 
rate when gas is draining on to the star and periods of low accretion rate when gas 
is accumulating in the inner disk. The instability could be a classical thermal insta¬ 
bility (50), of the type accepted as causing dwarf nova outbursts 111931 . or a related 
instability of dead zone structure 1 1120 .. 18]. The second category appeal to triggers 
independent of the inner disk to initiate the outburst. Ideas in this class are various 
and include perturbations from binary companions (59), the tidal disruption of radi¬ 
ally migrating gas clumps / giant planets 134211242) , and disk variability linked to 
a stellar magnetic cycle ini. Neither category of ideas is fully compelling (in the 
sense of being both fully worked out and consistent with currently accepted disk 
physics), so our discussion here will focus on a few key concepts that are useful for 
understanding current and future models of episodic accretion. 


5.1 Secular disk instabilities 

The classical instabilities that may afflict thin accretion disks are the thermal and 
viscous instabilities (276). These are quite distinct from the basic fluid dynamical 
instabilities (the MRI, the VSI etc) that we discussed in : (4] in that they address 
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the stability of derived disk models rather than the fluid per se. Thus the thermal 
instability is an instability of the equilibrium vertical structure of the disk, while 
viscous instability in an instability of an (assumed) smooth radial structure under 
viscous evolution. 

Before discussing how thermal or viscous instability might arise, we first define 
what these terms mean. Consider an annulus of the disk that is initially in hydrostatic 
and thermal equilibrium, such that the heating rate Q + matches the cooling rate Q . 
The heating rate per unit area depends upon the central temperature (equation [80]), 
and can be written assuming the a-prescription as, 

Q+= 9 vEQ 2 = 9 a k ^EQ. (179) 

4 4 jlniH 


The cooling rate directly depends upon the effective temperature, T^, but this can 
always be rewritten in terms of T c using a calculation of the vertical thermal structure 
(j p.3| l. In the simple case when the disk is optically thick, for example, we have from 
equation (87 1 that 7 f /7^ sk ~ (3t/4), and hence. 


e- = -r/. 

3t c 


(180) 


Both a and T may be functions of T c . Now consider perturbing the central tem¬ 
perature on a time scale that is long compared to the dynamical time scale (so that 
hydrostatic equilibrium holds) but short compared to the viscous time scale (so that 
E remains fixed j^] The disk will be unstable to runaway heating if an upward per¬ 
turbation to the temperature increases the heating rate more than it increases the 
cooling rate, i.e. if, 


dlog Q+ dlog Q 
dlogT c dlog T c ' 


The same criterion predicts runaway cooling in the event of a downward perturba¬ 
tion. A disk that is unstable in this sense is described as thermally unstable. It would 
heat up (or cool down) until it finds a new structure in which heating and cooling 
again balance. 

To determine the condition for viscous stability, we start by considering a steady- 
state solution E(r) to the diffusive disk evolution equation (52 1 . Following Pringle 
12761 we write ju = vE and consider perturbations ju —)• /i + 5/x on a time scale 
long enough that both hydrostatic and thermal equilibrium hold (in this limit T c is 
uniquely determined and v = V'(T)). Substituting in equation (52 1 the perturbation 
<5/1 obeys, 


<9 . <9/i 3 <9 


4/2 Z. 


(rV 2 8ii 


(182) 


The perturbation <5/i will grow if the diffusion coefficient, which is proportional to 
dfj./dE, is negative. This is viscous instability, and it occurs if. 


9 This may seem to require fine tuning, but in fact the ordering of time scales in a geometrically 
thin disk always allows for such a choice 12761 . 
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(183) 


A disk that is viscously unstable would tend to break up into rings, whose amplitude 
would presumably be limited by the onset of fluid instabilities that could be thought 
of as modifying the v(£) relation. 


5.1.1 The S-curve: a toy model 


The instabilities of interest for YSO episodic accretion can broadly be considered to 
be thermal-type instabilities. Noting that Q + aT c , and Q 7) 4 /r °c T t 4 /K (where 
K is the opacity at temperature T c ), we see that instability may occur according to 
equation ( | 18 1[ > if, 

• dlog Q + /dlog T c is large, i.e. if a is strongly increasing with temperature. 

• dlog Q /dlog T c is small, i.e. if K is strongly increasing with temperature. 

We expect a to increase rapidly with T c at temperatures around 10 3 K, as we transi¬ 
tion between damped non-ideal MHD turbulence at low temperature and the more 
vigorous ideal MHD MRI at higher temperature. We expect K to increase most 
strongly at temperatures around 10 4 K, as hydrogen is becoming ionized and there 
is a strong contribution to the opacity from H scattering (in this regime K can vary 
as something like r 10 Eol). Either of these changes can result in instability. 

Before detailing the specifics of possible thermal and dead zone instabilities in 
protoplanetary disks, it is useful to analyze a toy model that displays their essential 
features. We consider an optically thick disk, described by the usual classical equa¬ 
tions 11111 . whose angular momentum transport efficiency a and opacity K are both 
piece-wise constant functions of central temperature T c . Specifically, 


T c < Tcrit • 0! — Q'low i K — K"iow j 

T c A 7crit • tt = CZhigh, ^ = K*high> (184) 


with ai ow < rz hl „ h and fq ow < K'h; g ^. Our goal is to calculate the explicit form of the 
M(E) relation in the “low” and “high” states below and above the critical tempera¬ 
ture Tcrit- For a steady-state disk, at r 3> R*, heated entirely by viscous dissipation, 
the equations we need (mostly from jQ read. 
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a k B T c 

El jitw# 


(185) 


Note that the stellar mass M, and radius in the disk r enter these formulae only via 
their combination in the Keplerian angular velocity Q. Eliminating T c , 7 dis ^ t and 
V between these equations, we obtain a solution for M(E), 



valid on either the low or the high branch when the appropriate values for a and K 
are inserted. A solution on the low branch is possible provided that E < £ max , where 
A„iax is defined by the condition that T c = T cnl . Similarly, a high branch solution 
requires E > T mm with T c = 7 avt at T mm . The limiting surface densities are given 

by. 



(187) 


If K*high > Now and / or o: h ig h > ai ow , then Z max > E min and there will be a range of 
surface densities where accretion rates corresponding to either the low or the high 
branch are possible. 

Figure [23] shows, for a fairly arbitrary choice of the model parameters, the ther¬ 
mal equilibrium solutions that correspond to the low and high states of the disk 
annulus. One should not take the results of such a toy model very seriously, but it 
captures several features of more realistic models, 

• The solution has stable thermal equilibrium solutions on two branches, a low 
state branch where M for a given surface density is small, and a high state branch 
where it is substantially larger. In the toy model these branches are entirely sepa¬ 
rate, but in more complete models they are linked by an unstable middle branch 
(giving the plot the appearance of an “S”-curve). 

• There is a range of surface densities for which either solution is possible. 

• There is a band of accretion rates for which no stable equilibrium solutions exist. 

• The position of the S-curve in the E — M plane is a function of radius, with the 
band of forbidden accretion rates moving to higher M further from the star. 

The S-curve is derived from a local analysis, and the existence of annuli whose 
thermal equilibrium solutions have this morphology is a necessary but not sufficient 
condition for a global disk outburst. That time dependent behavior of some sort is 
inevitable can be seen by supposing that the annulus at 0.25 AU in Figure[23]is fed 
with gas from outside at a rate that falls into the forbidden band. No stable thermal 
equilibrium solution with this accretion rate is possible. If the disk is initially on the 
lower branch, the rate of gas supply exceeds the transport rate through the annulus. 
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Fig. 23 Example S-curves in the accretion rate-surface density plane from the toy model described 
in the text. For these curves we take Kj ow = Kj,igh = 1 cm 2 g~', ai ow = 10~ 4 , ahigh = 10~ 2 , and 
r crit = 10 3 K. The lower of the two curves is for = 1.6 x 10~ 6 s _1 (0.25 AU for a Solar-mass 
star), the upper for Q. — 5.6 x 10~ 7 s _1 (0.5 AU). 


and the surface density increases. This continues until E reaches and exceeds E max , 
at which point the only available solution lies on the high branch at much higher 
accretion rate. The annulus transitions to the high branch, where the transport rate 
is now larger than the supply rate, and the surface density starts to drop. The cycle 
is completed when E falls below T mm , triggering a transition back to the low state. 

The evolution of a disk that is potentially unstable (i.e. one that has some annuli 
with S-curve thermal equilibria) is critically dependent upon the radial flow of mass 
and heat, which are the key extra ingredients needed if an unstable disk is to “or¬ 
ganize” itself and produce a long-lived outburst. To see this, imagine a disk where 
annuli outside rj are already on the high branch of the S-curve, while those inside 
remain on the low branch. The strong radial gradient of T c implies a similarly rapid 
change in v, which leads to a large mass flux from the annuli that are already in 
outburst toward those that remain quiescent (equation [53|. The resulting increase 
in surface density, along with the heat that goes with it, can push the neighboring 
annulus on to the high branch, initiating a propagating “thermal front” that triggers 
a large scale transition of the disk into an outburst state. 

The quantitative modeling of global disk evolution including these thermal pro¬ 
cesses is well-developed within the classical a disk formalism 111931 . For a minimal 
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model, all that is needed is to supplement the disk evolution equation ( [52] ) with a 
model for the vertical structure (conceptually as described for the toy model above) 
and an equation for the evolution of the central temperature. This takes the form 


I551II951. 


dT c Q+-Q 
dt c p Z 



(188) 


Here c p is the specific heat capacity at constant pressure, and M is the gas constant. 
The first term on the right hand side describes the direct heating and cooling due to 
viscosity and radiative losses, while the second and third terms describe PdV work 
and the advective transport of heat associated with radial mass flows. In general 
there should be additional terms to represent the radial flow of heat due to radiative 
and / or turbulent diffusion (these effects are small in most thin disk situations, but 
become large when there is an abrupt change in T c at a thermal front). The treatment 
of these additional terms is somewhat inconsistent in published models, though they 
can significantly impact the character of derived disk outbursts 12541 . 

5.1.2 Classical thermal instability 

In physical rather than toy models for episodic accretion a and K are smooth rather 
than discontinuous functions of the temperature. A local instability, with a resulting 
S-curve, occurs if one or both of these functions changes sufficiently rapidly with 
T c (so that equation |l81| is satisfied). No simple condition specifies when a disk that 
has some locally unstable annuli will generate well-defined global outbursts, but 
loosely speaking outbursts occur provided that the branches of the S-curve (and the 
values of Z max and E mm ) are well separated. 

The classical cause of disk thermal instability is the rapid increase in opacity as¬ 
sociated with the ionization of hydrogen, at T ~ 10 4 K. Around this temperature K 
can rise as steeply as T 10 , and the disk will invariably satisfy at least the condition 
for local thermal instability. The evolution of disks subject to a hydrogen ioniza¬ 
tion thermal instability was first investigated as a model for dwarf novae (eruptive 
disk systems in which a white dwarf accretes from a low-mass companion star) 
Ml361 [234tt . and subsequently applied to low-mass X-ray binaries. Thermal insta¬ 
bility models provide a generally good match to observations of outbursts in these 
systems (which are of shorter duration that YSO outbursts, and correspondingly 
better characterized empirically), and are accepted as the probable physical cause. 
Good fits to data require models to include not only the large change in K that is 
the cause of classical thermal instability, but also changes in a between the low and 
high branches of the S-curve. Typically ai ow < 10 2 , whereas Cthigh ~ 0.1 111741 . 
MHD simulations that include radiation transport have shown that the S-curve de¬ 
rived from a disk models can be approximately recovered as a consequence of the 
MRI, and that the difference in stress between the quiescent and outburst states may 
be attributable to the development of vertical convection within the hot disk 1 1491 . 

By eye, the light curves of FU Orionis events look quite similar to scaled versions 
of dwarf novae outbursts, so the success of thermal instability models in the latter 
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sphere makes them a strong candidate for YSO accretion outbursts. The central tem¬ 
perature of some outbursting FUOr disks, moreover, almost certainly does exceed 
the 10 4 K needed to ionize hydrogen, making it inevitable that thermal instability 
physics will play some role in the phenomenon. Detailed thermal instability models 
of FU Orionis events were constructed by Bell & Lin Il50l . who combined a one¬ 
dimensional (in r ) treatment of the global evolution with detailed a model vertical 
structure calculations. They were able to find periodic solutions that describe “self- 
regulated” disk outbursts (i.e., requiring no external perturbation or trigger), with 
the disk alternating between quiescent periods (with M = 10 s — 10 7 M 0 yr 1 ) 
and outburst states (with M > 10 4 M 0 yr -1 ). These properties, and the inferred 
outburst duration of ~ 10 2 yr, are in as good an agreement with observations as 
could reasonably be expected given the simplicity of the model. 

The weakness of classical thermal instability models as an explanation for YSO 
accretion outbursts is that they require unnatural choices of the viscosity parame¬ 
ter a m. Thermal instability is indisputably tied to the hydrogen ionization tem¬ 
perature, which exceeds even the mid-plane temperatures customarily attained in 
protoplanetary disks. The required temperatures can be reached (if at all) only ex¬ 
tremely close to the star, and the radial region affected by instability extends to no 
more than 0.1-0.2 AU. The viscous time scale on these scales is short, so matching 
the century-long outbursts seen in FUOrs requires a very weak viscosity — Bell & 
Lin l50l adopt cthigh = 10 -3 . This is at least an order of magnitude lower than the 
expected efficiency of MRI transport under ideal MHD conditions l89l 1309 -1. More¬ 
over, in the specific case of FU Orionis itself, radiative transfer models suggest that 
the region of high accretion rate during the outburst extends out to 0.5-1 AU I356L 
substantially larger than would be expected in the thermal instability scenario. 


5.1.3 Instabilities of dead zones 


A dead zone in which the MRI is suppressed by Ohmic resistivity (( 4.4.3| l supports 
a related type of instability whose high and low states are distinguished primarily by 
different values of a, rather than by the thermal instability’s different values of K. 
The origin of instability is clear within Gammie’s original conception 11 191 of the 
dead zone, which has a simple two-layer structure. The surface of the disk, ionized 
by X-rayp 3 ) is MRI-active and supports accretion with a local a ~ 10 -2 . Below 
a critical column density Ohmic resistivity completely damps MRI-induced turbu¬ 
lence (according to equation |170[ ), and the disk is dead with a = 0. This structure 
can be bistable if the surface density exceeds that of the ionized surface layer. The 
low accretion rate state corresponds to a cool, externally heated disk with a dead 
zone; the high accretion rate state to a hot thermally ionized disk at the same E. 

Martin & Lubow 1 228112301 have shown that the local physics of Ohmic dead 
zone instability can be analyzed in an manner closely analogous to thermal instabil¬ 
ity. The sole difference lies in the reason why the lower state ceases to exist above a 


10 


Cosmic rays in the original model, though this is an unimportant distinction. 
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critical surface density (Z max in Figure |23)l. For thermal instability Z max is set by the 
onset of ionization at the disk mid-plane, and the attendant rise in opacity. A sim¬ 
ple dead zone, however, does not get hotter with increasing E, because the heating 
(either from irradiation, or viscous dissipation in the ionized surface layer) occurs 
at low optical depth at a rate that is independent of the surface density. It is then not 
obvious how even a very thick dead zone can be heated above 10 3 K, “ignited”, and 
induced to transition to the high state of the S-curve. 

One way to trigger a jump to the high state is to postulate that some source of 
turbulence other than the MRI is present to heat the disk. Gammie 111201 suggested 
that inefficient transport through the ionized surface layer would lead to the build 
of mass in the dead zone ona until Q « 1 and self-gravity sets in (see §4.2[ ). We 
can readily estimate the properties we might expect for an instability triggered in 
this manner by the onset of self-gravity at small radii. Suppose that at 1 AU the 
temperature in the quiescent (externally irradiated) disk state is 150 K. Then to 
reach Q = 1 requires gas to build up until E ~ 7 x 10 4 g cm 2 , at which point the 
mass interior is M ~ Kr 2 E ~ 0.025 M 0 . This is comparable to the amount of mass 
accreted per major FU Orionis-like outburst. In the outburst state the disk will be 
moderately thick (say h/r = 0.3 (49]), and the viscous time scale (l/a£2)(h/r)~ 2 
works out to be about 200 yr for a = 0.01, again similar to inferred FUOr time 
scales. At this crude level of estimate, it therefore seems possible that a self-gravity 
triggered dead zone instability could be consistent both with the inferred size of 
the outbursting region 113561 and with theoretical best guesses as to the strength of 
angular momentum transport in ideal MHD conditions. 

Time-dependent models for outbursts arising from a dead zone instability were 
computed by Armitage et al. m, and subsequently by several groups in both one¬ 
dimensional B359113571122811254] and two-dimensional models (358112411 . The more 
recent studies show that a self-gravity triggered instability of an Ohmic dead zone 
can give rise to outbursts whose properties are broadly consistent with those of 
observed FUOrs. Instability persists even if there is a small residual viscosity within 
the dead zone Il23l 12301. which could arise hydrodynamically in response to the 
“stirring” from the overlying turbulent surface layer B105I . 


It is established that Gammie’s dead-zone structure, as originally postulated 
111 19l . can be unstable to the development of a limit cycle in which outbursts 
and quiescent intervals alternate. It is much less clear whether a more realistic 
(we think) inner disk structure, affected by all three non-ideal terms (14.4.4 1 , 
evolves in a similar way. Ambipolar diffusion can damp turbulence in the low 
density surface layer, while the Hall effect (and winds) can lead to significant 
laminar stresses. It remains to be shown that our best theoretical models for 
disks on AU scales are locally unstable in the same way as a simple Ohmic 
dead zone, and to investigate the type of global outbursts that any instabilities 
might yield. The Hall effect and / or winds could in principle also lead to 
entirely different types of eruptive behavior. 
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5.2 Triggered accretion outbursts 

Accretion variability, including (perhaps) the large scale outbursts of FUOrs, can 
also be triggered by processes largely independent of the inner disk itself. Stel¬ 
lar activity cycles, binaries with small periastron distances, and tidal disruption of 
gaseous clumps or planets may all contribute. 


5.2.1 Stellar activity cycles 


As discussed in i 3.2.2| the inner disk is expected 1117811 and observed ||60l to be 
disrupted by the stellar magnetosphere. The complex dynamics of the interaction 
between the field — which may be misaligned to the stellar spin axis and have 
non-dipolar components — and the disk 119011 is likely the dominant cause of T 
Tauri variability on time scales comparable to the stellar rotation period (i.e. days 
to weeks). If the strength of the field also varies systematically due to the presence 
of activity cycles analogous to the Solar cycle, these could trigger longer time scale 
(years to decades) accretion variability. The simplest mechanism is modulation of 
the magnetospheric radius across corotation 179). When the field is strong and r m > 
r co the linkage between the stellar field lines and the disks adds angular momentum 
to the d isk, impe ding accretion in the same way as gravitational torques from a 
binary (j 3.2.3p* Gas then accumulates just outside the magnetospheric radius, and 
can subsequently be accreted in a burst when the field weakens. 

The viability of such magnetically “gated” accretion as an origin for large scale 
variability is limited by the short viscous time scale of the disk at r « r m , which 
makes it hard to accumulate large masses of gas if the stellar fields are, as expected, 
of no more than kG strength. Models CD suggest that significant decade-long vari¬ 
ability could be associated with protostellar activity cycles, but there is no clear 
path to generating FU Orionis outbursts. Activity cycles are more promising as an 
explanation for lower amplitude, periodic EXor outbursts ||88l . We note that there 
has been little if any work on the possible interactions between time-variable stellar 
magnetic fields and an inner disk that has its own net field. 


5.2.2 Binaries 

An eccentric binary with an AU-scale periastron distance could funnel gas into the 
inner disk, increasing the accretion rate and leading to an outburst if the increased 
surface density is high enough to trigger thermal or dead zone instability. This mech¬ 
anism was proposed by Bonnell & Bastien 159) . and has subsequently been studied 
with higher resolution simulations [ 270[[l08l . Although there are some differences 
between the predicted outbursts and those observed (this is almost inevitable, as the 
limited sample of FUOrs is already quite diverse), it is clear that close encounters 


11 In compact object accretion, this is described as the “propellor” regime of accretion fi53l - 
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from binary or cluster companions induce episodes of substantially enhanced accre¬ 
tion. The obvious prediction — that FUOrs ought to be found with observable binary 
companions or preferentially associated with higher-density star forming regions — 
is neither confirmed nor ruled out given the small sample of known objects. 


5.2.3 Clump tidal disruption 

A final possibility is that accretion outbursts could be triggered by the tidal disrup¬ 
tion of a bound object (a planet or gas cloud) that migrates too close to the star. The 
necessary condition for this to occur is given by the usual argument for the Roche 
limit. If we consider a planet with radius R p and mass M p , orbiting a star of mass 
M, at radius r, the differential (tidal) gravitational force between the center of the 
planet and its surface is. 


_ GM* GM, _ 2GM, n 

Oidal — t 7 D « — 5 Rp- 

r l ( r + R p )- r J 


(189) 


Equating the tidal force to the planet’s own self-gravity, F se if = GM p /R 2 p , we find 
that tidal forces will disrupt the planet at a radius r^ai given approximately by. 


r tidal 



(190) 


An equivalent condition is that tidal disruption occurs when the mean density of the 
planet p < M*/r 3 . 

It is difficult to tidally disrupt a mature giant planet. A Jupiter mass planet has a 
radius of R p ~ 1.5 Rj at an age of 1 Myr 112261 . and will not be disrupted outside the 
photospheric radius of a typical young star. (Though such planets, if present in the 
inner disk, could alter the course of thermal or dead zone instability [ 82 ] 82071 .) If 
tidal disruption is to be relevant to episodic accretion we require, first, that the outer 
disk is commonly gravitationally unstable to fragmentation, and, second, that the 
clumps that form migrate rapidly inward (in the Type 1 regime discussed by Kley in 
this volume) without contracting too rapidly. Numerical evidence supports the idea 
that clump migration can be rapid 1342 8 47 , 70 , 8 235 1 . though it is at best unclear 
whether contraction can be deferred sufficiently to deliver clumps that would be 
tidally disrupted on sub-AU scales nm Assuming that these pre-conditions are 
satisfied, however, Nayakshin & Lodato 11242 II studied the tidal mass loss from the 
close-in planets and its impact on the disk. They found that the tidal disruption of 
~ 20 Rj clumps, interior to 0 . 1 AU, led to accretion outbursts consistent with the 
basic properties of FUOrs. 

The primary theoretical doubts about tidal disruption as a source of outbursts 
concern the relative rates of inward migration and clump contraction, which are 
both hard to calculate at substantially better than order of magnitude level. Observa- 
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tionally, this process would produce outbursts in systems whose disks were young, 
massive, and probably still being fed by envelope infall. 

Predicting the evolution of the inner disk on year-to-century time scales from 
first principles is extremely hard, and none of the mechanisms for episodic 
accretion discussed in this section has been subject to strict enough scrutiny to 
merit dogmatic conclusions. Each comes with its own theoretical caveats, and 
no single mechanism seems well-suited to explain the entire range of outburst 
behavior seen in YSOs 12111, Episodic accretion may be caused by a significant 
modification of a known process, a combination of known processes, or a 
process yet to be identified. 


6 Single particle evolution 


The evolution of solid particles within disks differs from that of gas because solid 
bodies are unaffected by pressure gradients but do experience aerodynamic forces. 
We discuss here how these differences affect the motion of single particles orbiting 
within the gas disk, and how we can describe the evolution of a “fluid” made up of 
small solid particles interacting aerodynamically with the gas. Issues such as the rate 
and outcome of particle collisions, that are central to early stage planet formation, 
are treated in the accompanying chapter by Kley, and elsewhere fl5) . 

The key parameter describing the aerodynamic coupling between solid particles 
and gas is the stopping time. For a particle of mass m that is moving with velocity 
Av relative to the local gas, the stopping time is defined as, 


mA v 

|Edrag| 


(191) 


where F t j rag is the magnitude of the aerodynamic drag force that acts in the oppo¬ 
site direction to Av. Very frequently, what matters most is how the stopping time 
compares to the orbital time at the location of the particle. We therefore define a 
dimensionless stopping time by multiplying t s by the orbital frequency f2 K , 


% s = t s n K . 


(192) 


The dimensionless stopping time is also called the Stokes number. 

For our immediate purposes it largely suffices to describe aerodynamic effects 
in terms of the stopping time, but eventually you will want to translate the results 
into concrete predictions for how particles of various sizes and material properties 
behave. This requires specifying whose form depends on the size of the par¬ 
ticle relative to the mean free path of gas molecules, and (in the fluid regime) on 
the Reynolds number of the flow around the particle 13491 . If the particle radius s is 
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small compared to the mean free path A (5 < 9A/4) the particle experiences Epstein 
drag, with a drag force, 

Fdrag = -yWvth^V. (193) 

Here p is the density of the surrounding gas, and the thermal speed of the molecules, 


i'th = 


8 k B T 
itflniH 


(194) 


is roughly the same as the sound speed. Because Epstein drag is proportional to the 
velocity difference Av (rather than the more familiar square of the velocity differ¬ 
ence), the stopping time is a function of the particle properties that is independent 
of the velocity difference. For a spherical particle of material density p m , 

ts= P --. (195) 

P ''a, 

The mean free path in protoplanetary disks is of the order of cm (larger in the outer 
disk), so the Epstein regime is relevant for particles that range from dust to those of 
small macroscopic dimensions. Drag laws appropriate for larger bodies, which fall 
into the Stokes regime of drag, are given by Whipple ll349il . 


6.1 Radial drift 


The most important consequence of aerodynamic forces is the phenomenon of ra¬ 


dial drift. In (2.1.2 we showed that radial pressure gradients result in a gas orbital 
velocity that differs from the Keplerian value by 50-100 m s~ (using typical disk 
parameters at 1 AU). Most commonly, the gas is partially supported against gravity 
by the pressure gradient, and so rotates more slowly than the Keplerian value. Let us 
consider how this velocity differential affects the evolution of solids in various lim¬ 
its, initially ignoring turbulence and ignoring the feedback of aerodynamic forces 
on the gas, and then including these processes. 


6.1.1 Particle drift without feedback 

Large bodies (z s 2> 1) orbit at close to the Keplerian speed, and the effect of gas 
on their evolution can be considered as a simple “headwind” if the disk is sub- 
Keplerian. Suppose that the gas orbits at a speed vk — Av, with Av <C Vk- The drag 
force |/'d ra g| = mAv£ 2k/ T s does work at a rate. 


E — — 1-fBrag | vat. 


( 196 ) 
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that leads to a change in the orbital energy E = 
of the orbit. Noting that, 

g GM.m do 
2 a 2 dr 


GMpn/2a, where a is the radius 

(197) 


and equating the two expressions for E, we find that the orbit decays at a speed 
v r = da /At that is given by, 



(198) 


The radial drift of large bodies is inversely proportional to their Stokes number. 

The simple headwind argument fails for small particles with T v <C 1, which in¬ 
stead are forced to orbit at the gas speed by the strong aerodynamic coupling. The 
particles do not feel the pressure gradient, so their non-Keplerian orbital motion 
results in a net radial force, 


F r _(v k -Av) 2 GAf* _ 2 v k Av 

9 — ■ 

m a a A a 

Equating this to the drag force for radial motion at speed v r , |Td rag |/m = v r Qa/x s , 
we find that radial drift for small particles occurs at the terminal drift speed. 


V r = —2T s Av. 


( 200 ) 


This is the speed relative to the gas, so for a disk that is accreting there is an addi¬ 
tional component given by the gas’ radial velocity. 

Intermediate-sized particles orbit at some speed between that of the gas and that 
given by the Keplerian velocity. To derive the general rate of radial drift 1346113221 . 
we consider a gas disk whose orbital velocity is. 


1 /2 

P(|>,gas = V K (1 ) 


( 201 ) 


The parameter rj ( h/r ) 2 . For e xampl e, if the disk has X ^ r 1 and central temper¬ 
ature T c r^ 1 / 2 , we showed in (2.1.2 that r\ = (11/4 )(h/r) 2 . Defining the particle 
radial and azimuthal velocities to be vy and vy, respectively, the equations of motion 


are, 


diy 
"d T 




(t> v 0 igas ). 


( 202 ) 

(203) 


The azimuthal equation can be simplified by noting that the specific angular mo¬ 
mentum remains close to Keplerian, 


d 

dr 



— v r — (rv K ) = \v r v K - 
dr 2 


( 204 ) 
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Fig. 24 Particles drift inwards in a disk wherever dP/dr, due to the combination of aerodynamic 
forces and sub-Keplerian gas rotation. The radial drift time scale = r/|v r |, in units of the local 
orbital period, is plotted as a function of the dimensionless stopping time T y . The fastest radial drift 
occurs for z s = 1. The specific numbers shown in the plot are appropriate for a disk with h/r = 0.03 
and a values of 10~ 2 (lowest curve on the left-hand side of the figure), 10~ 3 and 10~ 4 . 


This yields. 


1 t s v r v K 

'V gas ~ 2 r 


(205) 


We now substitute for Q.^ in the radial equation using equation (2011. Discarding 
higher order terms we obtain. 


dv r 

dt 



(206) 


The dv r /dt term is negligible. Dropping that, we eliminate (v^ 
equations (|205[) and (206 [l to obtain. 


(r/\’K)t s Vgas -r/VK 
(vK/r)t s + (r/v K )t^ 1 


gas ) between 


(207) 


In terms of the Stokes number the final result for the particle radial velocity is, 


V r = 


T s ' ty.gas fj V K 

Ts + rr 1 


(208) 


The previously derived results for very small and very large particles are recovered 
by taking the appropriate limits. 
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The speed of the radial drift implied by equation ( |208[ ) is shown in Figure [24] 
Very rapid drift is predicted for particles with x s ~ 1. For our fiducial disk model 
with E r -1 , T c r -1 / 2 and h/r = 0.03, the radial drift time scale thrift = r/|v r ] is 
just 10 3 orbital periods — one thousand years at 1 AU! Particles with this stopping 
time are very roughly meter-sized, and their rapid drift is the origin of the “meter¬ 
sized barrier” that severely constrains models of planetesimal formation [741 ; 1601. 


6.1.2 Drift with diffusion 


If the disk is turbulent, small dust particles that are aerodynamically well-coupled to 
the gas will diffuse radially and vertically. (The same physics will apply to trace gas 
species, such as water or CO molecules.) Diffusion will tend to equalize the concen¬ 
tration of the dust or trace gas species relative to the dominant gaseous component 
in the disk. To derive an equation for the evolution of the trace species 1811 . we ini¬ 
tially ignore any radial drift due to aerodynamic effects and write the concentration 
of the trace gas or dust (generically the “contaminant”) as, 

/=Y’ (209) 

where E<j is the surface density of the contaminant. If the contaminant is neither 
created nor destroyed within the region of the disk under consideration, continuity 
demands that, 

+ V • = 0, (210) 

where F,./, the flux, can be decomposed into two parts: an advective piece describing 
transport of the dust or gas with the mean disk flow, and a diffusive piece describing 
the tendency of turbulence to equalize the concentration of the contaminant across 
the disk. For / < 1 we can assume that the diffusive properties of the disk depend 
only on the gas surface density, in which case the flux can be written as, 

F d =E d v-DZv(j±y (211) 


Here v is the mean velocity of gas in the disk and D is the turbulent diffusion co¬ 
efficient. The diffusive term vanishes if / is constant. Combining this equation with 
the continuity equation for the gaseous component, we obtain an evolution equation 
for / in an axisymmetric disk. In cylindrical polar co-ordinates. 


df 

dt 


J_d_ 
rE dr 




( 212 ) 


This result has the form of an advection-diffusion equation, with the advective com¬ 
ponent being due to the radial flow of the disk gas. It is easy to generalize this equa- 
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tion to account for the radial drift of larger particles that are imperfectly coupled to 
the gas, by adding an additional flux representing the radial drift speed El- 

Determining what is the appropriate value for the turbulent diffusion co¬ 
efficient D involves many of the same uncertainties that afflict the determi¬ 
nation of the turbulent viscosity. For trace gas species and very small dust 
particles the zeroth-order expectation is that D ss v 12391 . though simula¬ 
tions of non-ideal MHD disk turbulence show both significant deviations and 
anisotropy between the vertical and radial directions 13628 . For larger bod¬ 
ies there is a well-determined analytic scaling with the Stokes number of the 
particles 13531 . 


6.1.3 Particle pile-up 


Solids that experience significant radial drift tend to become concentrated (“pile- 
up”) in the inner disk ( 35411351 1. The basic effect is present in the simplest case 
where diffusion and inward drag by the mean flow are small, and feedback of the 
particles on the gas can be neglected. The radial drift speed is then, approximately. 


Vr- -fi t]v K , 


(213) 


with 1 ] oc ( h/r ) 2 . For particles that are in the Epstein regime of drag, the stopping 
time in the mid-plane is (from equations [195] and [5J, 


_ Pm S _ Pmh S 
PO V t h E V th ' 


(214) 


Since h = c s /Q., and c s and v,;, differ only by a numerical factor, we obtain. 


Ts 


K p m 

iY s 


(215) 


Suppose now (perhaps not very realistically) that the surface density profile of solids 
has attained a steady-state, such that the mass flux is constant with radius. Then, 


M c / = —I'KrE^Vr = constant, 


(216) 


and substituting for v r we find, 


Ed 

E 


OC 



(217) 


For a disk with constant (h/r) the steady-state concentration of solids increases 
closer to the star as r 1 / 2 . In the more realistic case of a flaring disk with mid-plane 
temperature profile, say, T c r l,/2 , the effect is stronger. A constant a model of 
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such a disk has a steady-state gas surface density profile E <x r 1 , with the solids 
following Ed « r~ 2 . 

6.1.4 The Nakagawa-Sekiya-Hayashi equilibrium 

Up till now we have implicitly assumed that the evolution of the gas is unaffected 
by the evolution of the solids within it. Obviously this can never be strictly correct. 
If a population of solid particles are losing angular momentum to the gas through 
aerodynamic forces and spiraling inward, the gas must gain a corresponding amount 
of angular momentum. If the solid to gas ratio is only of the order of 1%, however, 
one might suppose that the effect of the angular momentum exchange on the gas 
would be small and, perhaps, ignorable. This is only partially true. First, a number 
of processes, including vertical particle settling (97], gas loss via photo-evaporation 
1 3271141 or MHD winds, and radial drift itself l35ll . can boost the solid to gas ratio, 
at least locally. Second, the equilibrium solution for radial drift in the presence of 
back reaction on to the gas can be unstable to the streaming instability 0521 . which 
can result in strong localized clumping of the solids. 

The generalization of the radial drift formula (equation |208| l to account for the 
back reaction of the drift on the gas is known as the Nakagawa-Sekiya-Hayashi 
(NSH) equilibrium 12411 . The NSH solution is derived by considering the inter¬ 
action between solids and gas in a simple disk model that ignores the effects of 
turbulence and vertical gravity. Both the gas, with density p„, pressure P and veloc¬ 
ity v g , and the solid particles, with density density p p and velocity \ p , are treated as 
fluids that interact with each other via aerodynamic drag. They obey continuity and 
momentum equations of the form 13521 . 



(218) 

(219) 


( 220 ) 


( 221 ) 


We have replaced the continuity equation for the gas with the condition for incom¬ 
pressibility, which is valid provided that velocities remain highly subsonic. The only 
other differences between the equations for the two species are the presence of a 
pressure gradient term in the gas momentum equation, and the pre-factor in the 
aerodynamic drag term expressing the differing inertia of the two fluids. 

It is straightforward to derive the steady-state axisymmetric solution for the drift 
of solids and gas from the above equations 1 24111352 1. To illustrate the difference 
between the NSH and no-feedback solutions, we quote here just the result for the 
relative radial velocity of the solids and the gas. With our definition of t] (equa- 
tion|201|i this takes the form. 
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Fig. 25 The relative velocity of solids compared to gas (normalized by r\vg, the parameter specify¬ 
ing departure of the gas disk from Keplerian rotation) is plotted as a function of the Stokes number 
T t . From left to right the curves show the NSH equilibrium solutions for p p /p g = 10~ 2 (effectively 
identical to the no-feedback solution), 3 x 10 2 , 0.1, 0.3. 1, 3 and 10. 


W.rel — 


Pg 1] T Vg 

P 1 + (*sPg/p) 2 ' 


( 222 ) 


where p = p g + p p . Equation (2081 is recovered in the case where p p <C p g and 


feedback can be neglected (note that we have ignored any radial gas motions due to 
processes other than particle-gas coupling in this version of the NSH solution). 

The predicted relative radial velocity between the solid and gas components is 
plotted in Figure[25]for various values of the solid to gas ratio. For p p /p s - 1 0 2 — 
10 1 the rate of drift at a given value of the stopping time is very similar to the 
no-feedback solution given by equation (208 1 . For large values of p p /p g , however, 
the peak of the relative velocity curve (which in this limit is predominantly outward 
motion of the gas) shifts to somewhat higher values of z s . 


6.2 Vertical settling 

Aerodynamic drag also modifies the vertical distribution of solids relative to gas. 
We first consider the forces acting on a small particle of mass m at height z, above 
the mid-plane of a laminar disk. The vertical component of stellar gravity exerts a 
force. 


fgt'av — IllflfcZ 


( 223 ) 
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The gas in the disk is supported against this force by the vertical pressure gradient, 
but solid particles are not. If started at rest a particle will accelerate until the gravi¬ 
tational force is balanced by drag. In the Epstein regime (given by equation |l93| l the 
resulting terminal velocity is, 


t'settle — 


Pm 

P 



(224) 


Using numerical values appropriate for a 1 pm particle at z ~ h at 1 AU (p = 
6 x 1CU 10 g cm~ 3 , z = 3 x 10 11 cm, r t h = 10 5 cm s” 1 ) gives a settling speed v sett i e i=s 
0.06 cm s -1 . The settling time, defined as, 


z 

'settle — ] r 

| Vsettle | 


(225) 


is about 1.5 x 10 s yr. In the absence of turbulence micron-sized particles ought to 
settle out of the upper layers of the disk on a time scale that is shorter than the disk 
lifetime. 

Turbulent diffusion acts to counteract the effects of settling. If the particle fluid 
with density p p behaves as a trace species (i.e. p p /p -C 1) then it obeys an advection- 
diffusion equation l97l! 1151. 


dp P 

dt 





+ (&hsPpz) • 


(226) 


Steady-state solutions to this equation can be found in the limit where the particle 
layer is thin enough that the gas density is approximately constant across the particle 
scale height. The dimensionless friction time T s is then independent of z and we find, 


Pp 

P 


( Pp N 


z 2 ' 

VP j 

exp 

z=0 

2h2 p. 


(227) 


where h p , the scale height of the particle concentration p p /p, is, 


h p — 



(228) 


If the turbulent diffusivity is comparable to the turbulent viscosity, i.e. D ~ v, the 
ratio of the concentration scale height to the gas scale height is just. 


hp 

h 



(229) 


Solid particles become strongly concentrated toward the disk mid-plane whenever 
their dimensionless friction time substantially exceeds a. For reasonable values of 
a this requires substantial particle growth. 
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6.3 Streaming instability 

Youdin & Goodman 113521 demonstrated that the aerodynamically coupled system of 
gas and solids, described by the NSH equilibrium, is linearly unstable to the growth 
of perturbations. The instability, known as the streaming instability by (very rough) 
analogy with the two-stream instability of plasmas, taps the free energy present 
in the relative motion between the solid and gaseous fluids, which is ultimately 
sustained by the background gradient in the pressure. It provides the best-studied 
route to forming planetesimals (km or larger bodies that are largely decoupled from 
the gas) from smaller solids with T s ~ 1 or less. 

The pre-requisites for the existence of the streaming instability are two-way aero¬ 
dynamic coupling between gas and dust within a rotating system (with shear and 
Coriolis force). A minimal mathematical analysis can be performed in the “terminal 
velocity approximation”, in which the relative velocity between gas and dust is, 

VP 

Av =- t s . (230) 

P 

This approximation yields a third-order dispersion relation 1 352 iTHj , while a more 
complete analysis (still neglecting vertical stratification) results in a sixth-order sys¬ 
tem. The linear growth rates, plotted in Youdin & Goodman 13521 . are functions 
of T s and the ratio of solid to gas density, p c i/p s ■ Growth is typically substan¬ 
tially slower than dynamical, with the most unstable modes having scales <C h. For 
Ti ~ 10 2 and Pp/Pt, ~ 0.1, for example, the linear growth time scale is a few hun¬ 
dred orbits. 

A simple physical (as opposed to mathematical) explanation of the streaming 
instability is frustratingly elusive. (Analogies to “traffic jams”, or to the drag reduc¬ 
ing properties of pelotons in bicycle races, are more relevant to the strong clustering 
that results from the streaming instability than to its existence as a linear instability.) 
The reader distressed by this state of affairs may seek solace in papers by Chiang 
& Youdin (74), and by Jacquet et al. ED, who discuss the origin of instabilities in 
simplified or related physical systems. 

The relationship between the saturated state of the streaming instability (which 
is of greatest interest when the fluctuations in particle density are very strongly 
non-linear) and the linear growth phase is non-obvious, and requires numerical sim¬ 
ulations. Starting with the work of Johansen & Youdin Il63l , several authors have 
quantified the outcome of the instability in protoplanetary disks (for recent exam¬ 
ples, see e.g. 13(1 161 350 1). Simplifying greatly, the streaming instability depends 
on the local solid to gas ratio (or metallicity, with super-Solar metallicities Z > 10 2 
being favored 1 165 , 30)) and on the magnitude of the deviation from Keplerian ve¬ 
locity T]\’k/c s (with small values of this parameter promoting clumping |3 1 1 ). Using 
two-dimensional simulations Carrera et al. l67l suggest that the streaming instabil¬ 
ity can result in strong clumping for 10 3 < T s < 3, with the sweet spot where the 
lowest metallicity is required occurring for T s ~ 0.1 at Z « 0.015. 
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Fig. 26 Illustration of some of the processes that can lead to streaming instability in the aerody- 
namically coupled particle-gas system. (Simulations of the streaming instability and gravitational 
collapse by Jake Simon.) 


The existence of a strongly inhomogenous distribution of solid particles has im¬ 
portant implications for particle growth and planetesimal formation, irrespective of 
whether the over-densities are strong or very strong. There is particular interest, 
however, in determining whether the streaming instability can yield over-densities 
that exceed the Roche density (§5.2.3|>, given approximately by. 


Pp 


M* 


(231) 


Particle clumps whose density exceeds the Roche density can collapse gravitation¬ 
ally into planetesimals, whose properties (such as size and binarity 112441 ) will de¬ 
pend upon the statistics of the particle density field generated by the streaming in- 











Physical Processes in Protoplanetary Disks 


91 


stability. Simulations that include self-gravity show that collapse is a likely outcome 
in regions of the disk where the streaming instability is strong 11162111611 . 

The circumstances that lead to the formation of streaming-unstable regions 
within protoplanetary disks, along with the outcome of the instability when it oc¬ 
curs, are by no means definitively established. Figure 26 illustrates the flavor of 
theoretical models now under investigation ltl60l . which invoke the single particle 
processes discussed in this section as essential elements. Vertical settling and radial 
drift, operating on particles that have grown through collisions f58l to be imperfectly 
coupled to the gas, act to enhance the local metallicity toward the values where the 
streaming instability would operate. Settling and pile-up, however, may not always 
be sufficient, and the next section is devoted to processes that can generate structure 
and additional enhancement in the local metallicity within the disk. 


7 Structure formation in protoplanetary disks 

Up until now we have largely assumed that the gas and dust in protoplanetary disks 
follow axisymmetric distributions, with smooth (but very probably different ED) 
radial profiles. This is an approximation, which is known to fail spectacularly in 
some observed systems. We have already discussed the possibility that disks may 
be warped ( §3.4| i, and warps are observed in systems including HD 142527 1681 . 
Other observational indications of non-trivial disk structure include, 

• Classification of a significant fraction of protoplanetary disks as transitional 
disks cm based on evidence of inner cavities in (at least) the dust distribu¬ 
tion. 

• Radial structure in molecular emission linked to the presence of ice lines, for 
example of CO 12811 . 

• Multiple rings of emission seen in high resolution mm / sub-mm observations of 
HLTau El. 

• Pronounced non-axisymmetric (“horseshoe”-shaped) sub-mm emission in sys¬ 
tems including IRS 48 13391 . 

• Spiral arms and other non-axisymmetric structures seen in scattered light images 
of disks 11301 . 

An open and important question is whether these structures are a consequence of 
— or a precursor to — planet formation. That question cannot yet be answered, 
but keeping it in mind we discuss here a number of processes that can lead to the 
formation of directly observable (and hence necessarily large-scale) structure in one 
or both of the gas and dust distributions within disks. Independent of the topical 
observational interest, any process that can generate inhomogeneity in the solid dis¬ 
tribution is potentially important theoretically. Forming planetesimals appears to be 
hard (at least for us, though perhaps not in Nature), but could be made easier if there 
are processes that locally enhance the ratio of solids relative to the gas. 
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7.1 Ice lines 


The water snow line, together with the silicate sublimation front and various ice lines 
in the outer disk, are potentially critical locations for planet formation. Most often, 
this importance is quantified by noting that the equilibrium chemical composition of 
a Solar abundance gas has a substantially larger mass of condensible solids outside 
the snow line than inside (by about a factor of 4 in the classical Minimum Mass Solar 
Nebula fl42l . rather less than that using more modern calculations of the chemical 
equilibrium 82111 !. The likelihood that this leads to a jump in solid surface density at 
the snow line is then invoked as the reason why the Solar System has only terrestrial 
planets at smaller radii, and giants beyond. 

These arguments are valid but incomplete. First, the equilibrium chemical com¬ 
position is only linked directly to the solid surface density in the limit where the 
solid particles remain small and well-coupled to the gas. Particles that grow to be 
large enough that radial drift becomes significant will instead develop a surface 
density profile that is both different from that of the gas (16. 010531), and depen¬ 
dent on the size distribution. If icy particles are typically substantially larger than 
silicates (as is frequently suggested) their more rapid radial drift could lead to an 
instantaneously lower surface density of solids outside the snow line than inside. 
Second, although it is possible to construct models in which an assumed jump in 
planetesimal surface density at the snow line contributes to efficient core formation, 
the compositional effect is not the whole story. The greater area of planetary feed¬ 
ing zones at larger radii, along with more complex effects such as Type 1 migration 
and pebble accretion, affect the outcome of planet formation at different radii to a 
similar extent. In my opinion the most important role of ice lines may instead be as 
a preferential site for planetesimal formation, or perhaps as a location where Type 1 
migration stalls. 

The pressure in the protoplanetary disk is substantially below that of the triple 
point of water, and hence the snow line marks a radial transition between ice and 
water vapor. Under mid-plane conditions, the corresponding temperature is typi¬ 
cally T = 150 — 180 K. Where this isotherm lies in the disk is a function of the 
stellar luminosity and accretion rate — neither of which are constant over time — 
and of the disk opacity which may change due to coagulation. Theoretical models 
[ 1941 [12211236 1 suggest that r snow moves inward from « 3 AU, at such time as the 
accretion rate M « 10 7 M 0 yr 1 , to within 1 AU when the accretion has dropped 
to M ss 10 9 Mq yr _I . At still lower M the inner disk becomes optically thin, and 
the resultant rise in temperature pushes the snow line back out to 2-3 AU. 


Martin & Livio 82271 showed that the above estimates, calculated within the 
framework of relatively simple disk models that include viscous heating and 
irradiation, are significantly modified if the true disk structure instead resem¬ 
bles Gammie’s 81191 layered model with a mid-plane dead zone. The calcu¬ 
lation of snow line evolution may require further revision if winds — which 
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Fig. 27 Illustration of some of the physical processes occurring near ice lines in the protoplanetary 
disk. Icy materials drifting radially inward sublimate when they reach the ice line, releasing any 
higher temperature materials that were embedded into aggregates (76). The resulting vapor flows 
toward the star at the same speed as the rest of the gas, but also diffuses outward down the steep 
gradient in concentration (317) . It may then recondense, either into new particles or on to pre¬ 
existing particles (296) . Some combination of these effects may feedback upon the gas physics via 
changes to either the opacity or, in models where MHD processes dominate angular momentum 
transport, the ionization state HU). 


potentially change both the surface density profile and the fraction of potential 
energy that goes into disk heating — are important on AU-scales. 

Figure [27]illustrates some of the key physical processes occurring near ice lines. 
Where ice lines occur can be calculated by use of the Clausius-Clapeyron relation, 
which gives the saturated vapor pressure P e q at temperature T in terms of the latent 
heat L of the phase transition, 

P^=C L e ~ L ! m . (232) 

Here Si is the gas constant, and Q, is a constant that depends upon the species 
involved. For water, LjSi = 6062 K and Cl = 1.14 x 10 13 g cm -1 s“ 2 J76). We can 
compare this pressure to the actual partial pressure of vapor in the disk. Using water 
with molecular weight /J.h -,0 = 18 as an example, if the surface density of vapor is 
E v , the mid-plane pressure is, 
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Fig. 28 Example steady-state profiles of gas (upper solid line), radially drifting icy particles (solid 
blue lines) and water vapor (red dashed lines) in a turbulent protoplanetary disk. The assumed 
disk model has an accretion rate of 10~ 8 M a yr -1 , a temperature T = 150(r/3 AU) v 2 , and an 
a parameter of 5 x 10~ 3 . The ratio of the turbulent diffusivity to the turbulent viscosity is taken 
to be unity, and the concentration of icy solids is set (arbitrarily) to be 10~ 2 at 10 AU. The rapid 
radial drift of cm-sized particles leads to a high concentration of water vapor in the inner disk. 
Outward diffusion and re-condensation of this vapor — assumed here to form particles of a single 
fixed size — leads to an enhancement of solids just outside the snow line f3l7| . Note also the more 
elementary conclusion that the vapor concentration in the inner disk is a direct probe of the mass 
flux of radially drifting solids encountering the snow line. 


, = 1 gv k B T 

\p2n h Hn 2 oin H 


(233) 


If P v Pq q water ice will sublimate, whereas if / ( - Pq q vapor will condense into 
solid form. For small particles whose sizes are measured in mm or cm sublimation is 
rapid G6). and hence to a reasonable approximation sublimation and condensation 
processes act to maintain the vapor pressure close to the equilibrium value. 

The large value of L/2% implies that the snow line is a sharply defined transition 
within the disk. At 150 K, the characteristic temperature interval over which the 
equilibrium vapor pressure varies, P cq /(dP cq /dT), is just a few Kelvin. If sublima¬ 
tion of icy particles that drift through the snow line is fast, this has the consequence 
of imposing a sharp radial gradient in water vapor concentration at the snow line 
which, in a turbulent disk, will in turn result in a diffusive outward flux of vapor 
(equation |2 11 1 >. Re-condensation of the vapor can then lead to an enhancement of 
the solid surface density immediately outside the snow line. 

Stevenson & Lunine 13171 proposed that a vapor diffusion / condensation cycle 
of this type could enhance the surface density of ice enough to promote the rapid for- 
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mation of Jupiter’s core. The strength of the effect depends upon the size of the icy 
solids drifting in toward the snow line, and upon what is assumed about processes 
including disk turbulence, condensation and coagulation / fragmentation. Figure [28] 
shows the results of a particularly simple calculation (after l8Tll76 ;i), which assumes 
that vapor condenses (or, condenses and rapidly coagulates) into solid particles of 
a single size that matches the size of icy solids drifting in from larger radii. For 
the adopted disk model, mm-sized icy solids have relatively low drift velocities, and 
these particles sublimate into vapor without generating any local enhancement in the 
surface density of solids. Larger cm-sized particles, conversely, are enhanced by a 
factor of several outside the snow line as a consequence of the diffusive transport of 
vapor followed by condensation. Surface density enhancements of this magnitude 
could be important, particularly in models where planetesimal formation depends 
upon the disk locally exceeding a threshold value of metallicity (as is the case for 
the streaming instability, see Q (Ml). 

Ros & Johansen |296l investigated a related possibility. Instead of assuming that 
vapor condenses into new particles (or, on to very small grains released when ag¬ 
gregates break up at the snow line), they modeled the growth of pre-existing solids 
as vapor condenses on to their surfaces. A simple collisional argument gives the 
growth rate due to vapor condensation / sublimation as, 



(234) 


for a particle of mass m and radius 5, surrounded by vapor of density p,, and thermal 
speed i’ t h. (As noted by Supulver & Lin B319I . this is not an exact expression 111431 .) 
Using a Monte Carlo approach, Ros & Johansen 12961 found that condensation on 
to particle surfaces could provide an efficient growth mechanism up to sizes of the 
order of 10 cm. This could aid planetesimal formation by boosting the stopping time 
of particles into the range preferred by the streaming instability E). Moreover, by 
removing mass from the directly observable mm-size regime, condensation-driven 
growth could suppress the mm and sub-mm flux from disks in the vicinity of ice 
lines 13551 . 

The aforementioned physics affects only trace components of the disk — the icy 
solids and the resulting vapor. It is easy, however, to contemplate feedback processes 
that couple the evolution of solids at ice lines to the bulk of the gas disk. At a mini¬ 
mum the opacity will vary depending upon the radial distribution of solids. Beyond 
that, we have already noted that the efficiency of angular momentum transport (in 
MHD models) is expected to be a function of the local ionization state, and that 
the ionization balance is affected by the abundance of small grains. Kretke & Lin 
11851 suggested that the enhanced abundance of solids near the snow line would act 
to suppress the rate of angular momentum transport, and that this could lead to the 
formation of a local pressure maximum that would act to trap particles (producing, 
in principle, a positive feedback loop). This argument (which has been invoked in 
some models of collisional growth 1631 ) is highly plausible, though the breadth and 
complexity of the physics involved makes quantitative investigation challenging. 
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Although I have focused here on the snow line, analogous considerations carry 
over to other ice lines. The CO ice line corresponds to a temperature of T = 
17 — 19 K 05411 . and is somewhat more complicated to model because the CO is 
typically mixed with N 2 and water ices. The silicate “ice line” (or sublimation 
front) could also be important, since it lies at radii where a high fraction of 
stars are observed to host short-period planetary systems. 


7.2 Particle traps 


In a disk with a monotonically declining pressure profile, radial drift is always in¬ 
ward. More generally, however, aerodynamic drift is directed toward pressure max¬ 
ima, and can be outward if there is a local pressure maximum within the disk. This 
possibility was recognized in a prescient paper by Whipple 13491 . who appealed to 
it as part of a model for the formation of cornet^] The tendency for solids to be 
aerodynamically enhanced in the vicinity of pressure maxima is often described as 
particle “trapping”, though this term is somewhat misleading; in a turbulent disk 
small particles are at most temporarily detained by pressure maxima rather than 
being permanently trapped. 

The effect of local pressure maxima on the radial distribution of solids can be 
derived, in the limit where the solids remain a trace contaminant, following the 


methods described in 16.1.1 and 16.1.2 For an axisymmetric disk with an arbi¬ 
trary mid-plane pressure profile, the radial velocity v r of particles under the action 
of aerodynamic forces remains as given by equation ( |208[ ), with the parameter 77 
describing the deviation from Keplerian velocity becoming 13221 . 


(-V 

rdlnZ . 1 

I (n 


.1 1 v? J) 

din r 


(235) 


In this formula q is defined as the local power-law index describing the flaring of 
the disk, 

-« r q ~\ (236) 

r 

such that a non-flaring disk has q = 1. In axisymmetry and steady-state, equa¬ 
tions ( |210| i and (|2 1 1 [ > can be immediately integrated to give, 

K^diff + ^adv) =k, (237) 


where the diffusive and advective fluxes are. 


12 Quoting from his paper, “should it be possible for a toroid of higher density to occur in the Solar 
nebula, the growing planetesimals would be drawn toward it from the inside as well as from the 
outside...”. 

















Physical Processes in Protoplanetary Disks 


97 


-Fdiff = -DE 


d 

di- 



Tadv — V,'. 


(238) 


and the constant k is just the radial flux of solid material. Written out explicitly, the 
concentration of particles f=Ej/E obeys a first-order differential equation. 


d/_W k 1 

dr D J ~DEr' 


(239) 


Analytic solutions to this equation are possible for simple choices of v r , D and E. 
A straightforward quadrature gives the solution for the concentration profile given 
more realistic choices of these functions. 

Figure 29 shows an illustrative numerical solution to equation (239i. For this 
example, we have modeled a trap by assuming (arbitrarily) that the viscosity in the 
gas disk is reduced across a moderately narrow annulus. The lower viscosity leads to 
a higher surface density, producing a pressure maximum which in turn concentrates 
particles. The concentration effect is strongly size-dependent. Small particles (in this 
example those with radii of 0.1 mm and 1 mm) have a radial velocity that is similar 
to that of the gas, and are largely unaffected by the pressure maximum. Larger cm- 
sized particles, on the other hand, have a larger magnitude of radial drift, and can be 
strongly concentrated at the location of the pressure maximum. The enhancement in 
the local particle density can reach several orders of magnitude, depending both on 
the particle size and on the “strength” (radial width and amplitude) of the pressure 
maximum within the gas disk. 


As should be obvious from the preceding discussion, in a turbulent disk (and 
ignoring particle feedback on the gas) it is always possible to find a steady- 
state solution for the radial particle concentration in which the particle mass 
flux is a constant at all radii. Physically, this is because particles accumu¬ 
late near pressure maxima until the radial gradient in concentration becomes 
large enough for turbulent diffusion to allow them to leak out 1 343 3601 , Pres¬ 
sure maxima only act as a “filter” 112901 , permanently removing large particles 
from the inward radial drift flow, when additional physical effects are included 
ll3601 . If the local concentration of large particles becomes large enough, for 
example, planetesimal formation 111371 could cause some fraction of the solid 
material to drop out of the radial flow. 


A number of physical effects have been identified that could lead to the formation 
of local pressure maxima within protoplanetary disks. These include the exterior 
edges of planet-carved gaps ( 23911290113601 , the outer edges of cavities created by 
photoevaporation ( §8.1| |(4]), and a photoelectric heating instability that may operate 
in gas-poor systems (primarily debris disks) I53ll222 ’l. Spiral arms in self-gravitating 
disks 1 29211123 1 and the inner edges of dead zones 12211 can also concentrate solids 
via closely related physical processes, though these environments typically involve 
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Fig. 29 The steady-state radial distribution of solids in a turbulent disk with an axisymmetric local 
pressure maximum (a particle “trap”). Upper panel: the surface density of the gas. Middle panel: 
the radial velocity of 0.1 mm (red), 1 mm (green) and 1 cm (blue) particles. The smallest particles 
have a radial velocity that is almost indistinguishable from that of the gas, while the larger particles 
experience rapid radial drift that can be outward near the location of the pressure maximum. Lower 
panel: the concentration C = Td/Tgas- normalized to an arbitrary value of 10~ 2 at 100 AU. The 
assumed disk model for this calculation has M = 10~ 8 M 0 yr _1 , M, = M 0 , h/r = 0.5, a = 10~ 3 
and D = v (c.f. equation |2 1 2[ >. The trap is modeled as a gaussian-shaped reduction in a to a 
minimum of 10~ 4 , with a width of 4h. The particles are assumed to be spherical, with a material 
density of 1 g cm~ 3 , and to follow the Epstein drag law. 


significant non-axisymmetry. With the exception of self-gravity, these possible lo¬ 
cations for pressure maxima either form at specific places within the disk (e.g. at 
the inner dead zone edge, defined by a characteristic mid-plane temperature), or at 
a time after when we expect planets to have formed (during photoevaporative disk 
clearing, or in response to pre-existing massive planets). It is possible, however, for 
the turbulence within the gas disk to be generically unstable to the formation of 
pressure maxima within zonal flows. This would be interesting because it would 
imply the (possibly transient) existence of multiple particle traps within the disk, 
that could play a role in early-stage planet formation 1271 1 . 
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7.3 Zonal flows 


It is evident from equation ( fl2| that an equilibrium can be set up in which radial 
forces from a complex pressure profile (that may include local pressure maxima) 
are balanced by radial variations in azimuthal velocity. In a local description (equa¬ 
tion 127| but here ignoring magnetic fields) the balance is between the pressure 
gradient term VP/p and that describing the Coriolis force 2Qq x v. When these 
terms balance, such that 

^ = 2Q 0 pv y , (240) 

the system is said to be in geostrophic balance (here x is the radial direction, and 
y the azimuthal). The pressure gradient is compensated by variations in the orbital 
velocity, creating zonal flows analogous to the banded structure of winds in giant 
planet atmospheres. 

A disk zonal flow is an equilibrium solution to the fluid equations, but that equi¬ 
librium may not be stable. Too pronounced a deviation from Keplerian rotation re¬ 
sults in a shear profile that is unstable to Rossby wave instability ll212l . We will dis¬ 
cuss this instability, which is similar to Kelvin-Helmholtz instability, in §7.5| Even if 
the rotation profile is stable, the diffusive nature of a classical viscosity (equation[5l) 
would tend to erase any small-scale perturbations in the pressure that are sourced 
from surface density fluctuations. Persistent zonal flows are thus not expected in 
classical disk theory. They have been observed, however, in local numerical sim¬ 
ulations of MHD turbulent disks 1116411 . The key to their formation appears to be 
the ability of MHD disk turbulence to generate large-scale structure in the magnetic 
field 0091 . which could be viewed as an inverse cascade of turbulent power. The 
details of how and when zonal flows form are not entirely clear, though Johansen 
et al. ll 1641 describe a simplified dynamical model in which large-scale variations 
in the Maxwell stress lead first to azimuthal velocity perturbations and thence to 
axisymmetric structure in the pressure and density. 

The lifetime and radial scale of zonal flows in protoplanetary disks depend upon 
the same factors that determine the properties of MHD disk turbulence more gener¬ 


ally (14.4.4i, namely the strength of non-ideal terms in the induction equation and 
the presence of net vertical magnetic field. In ideal MHD, lifetimes of tens of orbital 
periods and radial scales of the order of 10 /; appear to be typical 1164. 309 1, though 
these results require a double dose of caveats — first because the inferred scales are 
not much smaller than the size of the local simulation domains used, and second 
because they are large enough that curvature terms neglected in local models may 
be important. Nonetheless, the amplitude, scale and lifetime of zonal flows under 
ideal MHD conditions plausibly lead to strong local enhancements in the dust to 
gas ratio for particles with stopping times t s ~ 0.1 — 1 li92l . The presence of net 
vertical fields substantially enhances the amplitude of zonal flows ll34ll . 

Zonal flows are also found in MHD disk simulations that include ambipolar dif¬ 
fusion, and can be comparable in amplitude to the ideal MHD case if the net field 
is sufficient to stimulate a significant a ~ 10 2 0061 . However, both inferences 
from local simulations 13061 , and explicit tracking of particles in global simulations 
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(3621, suggest that zonal flows in the outer regions of protoplanetary disks, where 
ambipolar diffusion dominates, have properties close to the boundary beyond which 
strong particle enhancement would be expected. 

In summary, zonal flows are likely to be present in the inner disk, where 
ideal MHD is a good approximation, though these flows would only strongly 
influence the dynamics of relatively large solid particles. In the outer disk, 
where ambipolar diffusion is important and even mm-sized particles have sig¬ 
nificant stopping times, zonal flows could introduce observable large-scale 
axisymmetric structure and may contribute to particle concentration. In the 
Hall-dominates regime that prevails around 1 AU theoretical expectations are 
less clear. Extremely strong zonal structures were observed in local vertically 
unstratified Hall-MHD simulations 11891 . whereas comparable stratified runs 
instead led to large-scale Maxwell stress jT95l . It is therefore unclear whether 
there are circumstances in which zonal flows on AU-scales could contribute 
to particle concentration and planetesimal formation. 


7.4 Vortices 


Few issues in planet formation are as long-debated as the possible role of vortices. 
Very general arguments suggest that large-scale vortices could be present in pro¬ 
toplanetary disks and play an important role in planetesimal formation. We note, 
first, that disks are (approximately) two dimensional fluid systems, that in contrast 
to three dimensional systems support an inverse cascade of turbulent energy toward 
large scales II181 II . Second, for a barotropic disk (i.e. P = P{p) only) the vortensity 
co/p is conserved (equation 113 i. Taken together, these properties imply that vor¬ 


tices within disks have the potential to form persistent long-lived structures. Indeed, 
simulations of strictly two dimensional flows show that disks seeded with small- 
scale vorticity perturbations evolve to form a small number of large and persistent 
anticyclonic vortices II1241II671 . Anticyclonic vortices are high pressure regions 
that attract marginally coupled solids l44l 13251 . potentially catalyzing the subse¬ 
quent formation of planetesimals. Even absent planetesimal formation, the natural 
tendency of vortices to form large-scale non-axisymmetric dust features makes it 
tempting to identify them with observed disk asymmetries 13391 . 

The basic properties of disk vortices are well-established. What is much trickier 
is to determine (1) whether vortices form spontaneously in disks or only after planet 
formation (for the reasons already mentioned, spontaneous vortex formation gener¬ 
ally requires non-barotropic processes), and (2) whether three dimensional instabil¬ 
ities and / or particle feedback are fatal impediments to their survival. Observations 
as well as theory are probably needed to resolve these issues. 
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7.4.1 The Kida solution 

The magnitude of the vorticity in a strictly Keplerian disk is C0k = —(3/2)12 k- A 
vortex can be modeled as a spatially localized elliptical perturbation, within which 
the vorticity ft) = g>k + ft) v , with (o v a constant. Other types of vortex are possible, but 
rather remarkably this type can be described by an exact non-linear solution 11731 
that is useful for both analytic and numerical studies. 

The Kida solution 11731 describes a vortex within a shearing-sheet approximation 
to disk flow. Following Lesur & Papaloizou 11981 we define a cartesian co-ordinate 
system (x,y) that is centered at radius ro and which co-rotates with the background 
disk flow at angular velocity Qk = (ro). 


^ = r 0 </> 
y = ~(r-ro). 


(241) 


Kida considered time-dependent vortex solutions, but here we will worry only about 
vortices that are stead}] 13 [ Time-independent solutions are possible if the semi-major 
axis of the vortex is aligned with the azimuthal direction (x in the shearing sheet 
model), and the vorticity perturbation satisfies, 



(242) 


Here x = o/b is the aspect ratio of the vortex, which forms an elliptical patch with 
semi-major axis a and semi-minor axis b. The right-hand-side is evidently positive, 
which implies that the only steady Kida vortices in Keplerian disks are anticyclonic 
(with ft),, having the opposite sign to fix)- 

The complete Kida solution is written in terms of a streamfunction t// in an ellip¬ 
tic co-ordinate systems (fi, v), where. 


x = f cosh(^)cos(v), 
y = /sinh(/i)sin(v), 


(243) 


and / = a\J(x 2 — l)/x 2 - The solution can be split into a core and an exterior part. 



+ ^—|exp[-2(^-^ 0 )]cos(2v) , (244) 


13 For a derivation of the steady Kida solution, see e.g. the appendix of Chavanis EU. 
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Fig. 30 Contours of the Kida vortex streamfunction y(x,y) are shown for different values of the 
vortex aspect ratio x=a/b. Within the vortex core, delineated by the bold contour, the streamlines 
defined by the solution are elliptical, with fixed aspect ratio. Outside the core, the vortex merges 
smoothly into the background shear flow of the disk. 



which match at /i = /io = tanh ] (x 1 )• The cartesian velocity field is then given by 
v x = —d\jf/dy, v y = d \\! jdx. In general the cartesian representation of the velocity 
has no simple form, but within the core it is, 


Fr.core — 
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(245) 


describing simple elliptical motion. Figure [30] shows the contours (logarithmically 
spaced) of the full streamfunction for Kida vortices of varying aspect ratio. 

Barge & Sommeria l44l studied the trajectories of aerodynamically coupled 
solids that encounter vortices (using a different and approximate vortex model). 
The high pressure within anticyclonic vortex cores acts as an attractor for solids. 
The cross-section for capture is maximized for particles with dimensionless stop¬ 
ping time Tv = 1 ED. Since vortices can potentially grow to have radial extents 
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Armh (the supersonic velocity perturbations of larger vortices would radiate sound 
waves sapping their energy) a single large vortex can potentially trap a substantial 
mass of solids flowing radially toward it due to ordinary radial drift. 


7.4.2 Stability of vortices 

Although the geometry of protoplanetary disks is approximately two dimensional, 
the fact that disk vortices are limited in radial size to Arm h means that they are three 
dimensional objects that notice the vertical stratification. Barranco & Marcus l46l 
and Shen et al. 0041 . using three dimensional simulations, found that mid-plane vor¬ 
tices that would be highly stable in 2D are rapidly destroyed by three dimensional 
instabilities. The origin of these observed instabilities, at least in part, appears to be 
the elliptical instability 117211 . which occurs whenever there is a resonance between 
the vortex rotation period and inertial waves within the disk. In a disk environment, 
Lesur & Papaloizou 111981 find that purely gaseous vortices are unstable for almost 
all choices of the vortex aspect ratio and degree of vertical stratification, though 
these parameters strongly affect the linear growth time of the instability. The in¬ 
stability that afflicts 3D vortices, however, is typically slow-growing and of small 
radial scale. Numerically this makes studies of vortex survival particularly challeng¬ 
ing. Physically it means that the questions of vortex formation and vortex survival 
are closely linked, except perhaps at very large radii (where “primordial vortices” 
might persist for interesting periods of time) the vortex population in a disk is ex¬ 
pected to reflect an equilibrium between formation and destruction processes. 

A significant loading of solids will also impact vortex longevity, generally for 
the worse. Railton & Papaloizou 12861 studied the stability of generalized Kida vor¬ 
tices, containing both gas and dust in the limit of strongly aerodynamically coupled 
particles. They found that these configurations were vulnerable to parametric insta¬ 
bilities in the same way as gas-only vortices. This is consistent with a wide range 
of other analytic and numerical work ll72l 1154 . 1171 , which suggests that dust to 
gas ratios in the range between 0.1 and 1 are sufficient to imperil the survival of 
disk vortices. This limitation may not, however, preclude vortices playing a role in 
planet formation. If we assume that there is a continual source of vorticity within the 
disk, then vortex particle concentration up to p p ~ p„ could be sufficient to initiate 
planetesimal formation via a variation of the streaming instability 12821 . 

Observations have identified a number of systems (primarily transition disks) 
that show a non-axisymmetric distribution of sub-mm emission |69 157 339 1, 
consistent with that expected if aerodynamically coupled solids are accumu¬ 
lating in a vortex 13611 . The putative vortices in these examples may all be 
caused by planets. It would be interesting to ask whether useful constraints on 
vortices could be derived from observations of the most-axisymmetric disks, 
where there is no independent suspicion that planets already exist. 
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7.5 Rossby wave instability 

The Rossby wave instability (RWI) li2T2[ 1202; 201 1 is a well-characterized mecha¬ 
nism for producing vortices within protoplanetary disks. The RWI is a linear insta¬ 
bility that grows whenever there is “sufficiently sharp” radial structure in the disk. 
Specifically, for a two dimensional disk model with angular velocity Q. (r), verti¬ 
cally integrated pressure P(r), and adiabatic index y, we define the entropy S and 
epicyclic frequency K via. 


£ 7 ’ 


(246) 


In terms of these quantities, the stability of the disk to the RWI is determined by the 
radial profile of a generalized potential vorticity. 


*(r) = 


2QE 


: S~ 2 ^. 


(247) 


A necessary condition for RWI is that have an extremum. A precise sufficient 
condition is not known, but variations in the potential vorticity of the order of 10% 
over radial scales k, h appear to be enough to trigger instability, which leads to 
the formation of anticyclonic vortices on time scales that can be rapid — of the 
order of 0.1X2 1 12021 . The instability, which can be understood in terms of the 
local trapping of waves in the vicinity of the vortensity perturbation [125) 3371 1, has 
similarities to the Papaloizou-Pringle 112641 instability of accretion tori. The RWI is 
essentially a two dimensional instability B 2321l205l . though the vortices that it forms 
are as vulnerable to unrelated three dimensional instabilities as any others. 

We have already remarked that the edges of dead zones (and ice lines 1118510 are 
places where local pressure maxima may form. The RWI criterion (equation |247[ ) is 
not necessarily satisfied at every local pressure maximum, but it is nonetheless true 
that dead zone edges are plausible locations where the RWI may occur. Hydrody¬ 
namic models 1 340 154 2201 . and MHD simulations that include Ohmic diffusion 
as a simple model of dead zones 1 223 , 103, 1224 . f~06l , support this expectation, and 
show that the radial disk structure introduced by dead zones is likely to be unstable 
to the RWI and subsequent formation of vortices. The vortices in turn act as sites 
of particle concentration H220L Depending upon the radii involved, a single vortex 
generated by the RWI may have a lifetime that is very short as compared to the disk 
lifetime (this will be particularly true at the inner edge of the dead zone). However, 
in this scenario fresh generations of vortices may be expected to form as long as the 
non-ideal disk physics that sustains the RWI-unstable dead zone structure persists. 

A second location where RWI may occur is at the edge of a gap created dynam¬ 
ically by a massive planet 111791 :901. Figure 31 shows the evolution of a disk in 
a two-dimensional simulation of this scenario. The formation of a massive planet 
creates an approximately axisymmetric gap within the disk, whose edges can be 
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Fig. 31 Snapshots showing the hydrodynamic evolution of an almost inviscid disk containing a 
massive planet EZ). The planet rapidly clears an annular gap within the disk, whose edge is unsta¬ 
ble to the generation of vortices. The system then evolves through a phase when the outer gap edge 
hosts a single large vortex, which can be an efficient trap for solid particles. The disk has a = 10~ 4 
and h/r = 0.05 at the location of the planet, which has a mass ratio to the star of 5 x 10~ 3 . 


unstable to the RWI. The vortices that are generated fairly rapidly merge, creating 
a single large vortex at either edge of the gap then can subsequently trap particles. 
This process works best if the disk in the vicinity of the planet’s orbit has a low 
viscosity (roughly a ~ 10 3 or lower), and is at least partially a transient effect — 
the vortices form during the phase when the planet accretes its gaseous envelope. 
For these reasons, observable structure from planet-initiated vortices is likely to be 
easiest to see in the outer disk, where ambipolar diffusion damps turbulence 13081 
and the absolute lifetime of a single generation of vortices against disruptive insta¬ 
bilities is long. The particles trapped efficiently in the outer disk include those with 
s ~ mm that can be seen in sub-mm observations of protoplanetary disks 1 136311361 1. 


8 Disk dispersal 

There is not obviously any leading order mystery as to where the gas in protoplane¬ 
tary disks goes to. If we divide the mean disk mass estimated from sub-mm studies 
in Taurus (Mj ~5x 10“ 3 M 0 (8j) by the median accretion rate (M ~ 10~ 8 M 0 yr _I 
11351 ) estimated in the same region, we obtain a characteristic evolution timescale 
of 0.5 Myr. The lifetime of detectable disks might be expected to be a small multiple 
of this timescale — say a few Myr — which is indeed what is observed 11471 . 

The above exercise establishes that, from a purely observational perspective, 
most or all of the gas in protoplanetary disks could be lost via accretion on to the star. 
Beyond that, it proves nothing. Disk mass estimates and measurements of stellar ac- 
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cretion rates are subject to uncertainties, which when combined probably allow us 
to shift the inferred characteristic lifetime by an order of magnitude in either direc¬ 
tion. In particular, if disk masses are systematically under-estimated from sub-mm 
continuum observations (because of particle growth to sizes too large to contribute 
to the sub-mm opacity) the characteristic lifetime could approach or exceed the ob¬ 
served one, challenging the view that most of the gas is accreted. 

The next order observational diagnostic probes the time-dependence of the de¬ 
cay of disk signatures. Here a puzzle does emerge. The classical disk evolution 
equation ( [52] ) admits a self-similar solution 112 191 that ought to approximate the 
evolution at sufficiently late times. For a disk with a viscosity that scales with radius 
as v oc r f , the late-time evolution of the surface density close to the star is predicted 
to follow (e.g. Ifl5l ), 

jocfP/M/P-rt (248) 

If y = 1, for example, as would be the case for a disk with a steady-state surface 
density profile E oc r 1 , then the predicted late-time decay goes as f~ 3 / 2 . This is 
relatively slow, and would probably lead to a population of stars with weak disk sig¬ 
natures (in gas and dust tracers) that are not observed (for a review of these obser¬ 
vational arguments, see Alexander et al. 0). One might argue, of course, that this 
is an illusory problem that we have created for ourselves by placing unwarranted 
faith in the validity of classical viscous disk theory. A resolution along these lines 
is possible. More commonly, however, the discrepancy between the simple theory 
and observations is taken to imply that some distinct process acts to rapidly disperse 
the disk and terminate accretion. Photoevaporation is the leading candidate, on ac¬ 
count both of robust theoretical estimates that suggest it should be important, and 
observations that are consistent with photoevaporation occurring in some relatively 
extreme situations. 


8.1 Photoevaporation 

Disk photoevaporation is a purely hydrodynamic process that occurs when molecu¬ 
lar gas in the disk is dissociated or ionized by high energy (UV or X-ray) photons. If 
the gas is heated sufficiently to become unbound, it accelerates away from the disk 
under the influence of pressure gradient forces to form a thermally driven wind. 
Early models for photoevaporative flows were developed in the 1980s by Bally & 
Scoville 1431 in the context of massive stars surrounded by neutral disks, and in more 
quantitative detail by Begelman et al. l48l who studied X-ray heated disks around 
compact objects. The essential physics is thus very well-established. We will begin 
by considering how photoevaporation works if the disk is exposed to extreme ultra¬ 
violet (EUV) photons that have sufficient energy (Irv > 13.6 eV) to ionize hydrogen. 
This is probably not the dominant driver of photoevaporation from protoplanetary 
disks around low mass stars, but it is amenable to an analytic treatment that exposes 
the main principles 1 1511 1. 
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8.1.1 Thermal winds from disks 


We consider a disk whose surfaces are illuminated by a source of high energy pho¬ 
tons, that may come either from the central star or from other luminous stars within 
a cluster. The radiation heats a surface layer of the disk to a temperature T, with a 
sound speed c s . A characteristic scale r g (the “gravitational radius”) can be defined 
by asking where the sound speed equal the local Keplerian velocity. 


CM, 


(249) 


Noting that the thermal energy per particle is ~ k/fk, the radius r„ is approximately 
equivalent to the smallest radius where the total energy of the gas in the heated 
surface layer (i.e. thermal plus gravitational) is zero. For radii r > r„ the total energy 
is positive, and it is energetically possible for the surface gas to flow away in a 
thermal wind. 

This is all quite rough, and we should really consider both the hydrodynamic 
structure of the wind (which is similar to the textbook example of a Parker wind 
13451 ) and the role of rotation 112041 . Doing so results in an improved estimate of 
the critical radius beyond which a thermal wind is launched, which scales with but 
is significantly smaller than r g . 


i 0.2 


GM* 


1.8 

\M e 


)( 


lOkms- 


-2 


AU. 


(250) 


We have picked 10 km s _1 as a fiducial sound speed because this is approximately 
the sound speed in EUV-ionized gas, which has a temperature near 10 4 K. X-rays or 
far ultraviolet (FUV) photons (those with 6 eV < hv < 13.6 eV that dissociate but 
do not ionize hydrogen) heat the surface to lower temperatures — more like 10 3 K 
— resulting in correspondingly larger critical radii. 

Interior to r g the hot gas is bound, and unless some other process intervenes (such 
as a stellar wind or an MHD disk wind) it will form a static is otherm al atmosphere 
with a scale height that varies with radius as h r 3 / 2 (see (2.1.1 1 . Outside r g it 


will flow away, at a speed of the order of the sound speed. In the case of EUV 
illumination the hot and ionized surface layer is separated from the underlying cool 
gas by a sharp ionization front, which gives a clearly defined “base” to the wind. If 
the number density at the base is no, we expect a mass loss rate per unit area of the 
disk that is given by, 

t w ~ 2nm H n 0 c s , (251) 

up to factors of order unity that depend again on the detailed hydrodynamic structure 
of the flow. In the EUV case the mass loss profile due to photoevaporation is then 
determined by the radial scaling of the base density iio(r). Noting that the integrated 
mass loss rate is just. 


M,., = 


r r out 

L : 


27trE,,.dr, 


(252) 
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Fig. 32 Illustration of the simplest model of internal photoevaporation driven by extreme ultravi¬ 
olet (EUV) radiation (based on the “weak stellar wind” case from Hollenbach et al. (l5ll ). Stellar 
EUV radiation ionizes and heats the surface layers of the disk. Where the thermal energy of the 
surface layer remains small compared to the binding energy, the hot gas forms a bound atmosphere. 
At larger radii, where the gas is more weakly bound, the hot gas flows away in a thermally driven 
wind. Details of the radiative transfer and heating processes differ depending upon the nature of 
the high energy radiation, but a qualitatively similar scenario applies also to X-ray and FUV-driven 
photoevaporation. 


we conclude that if no(r) declines more steeply that r 2 the mass loss is predomi¬ 
nantly from the inner disk (near r c ), whereas a shallower profile leads to most of the 
mass being lost from the outer disk. 

Actually finding no(r) for photoevaporation driven by EUV radiation from the 
central star requires the solution of a radiative transfer problem, whose geometry 
is illustrated in Figure [32] The base density is determined by equating the rate of 
ionization to that of recombination, which occurs at a rate per unit volume an e n p 
with a , the recombination co-efficient, given by a w 3 x 1CV 13 . The main diffi¬ 
culty is determining the radial scaling of the ionization rate. This in principle has 
two components, a “direct” flux from the star and a “diffuse” field that originates 
from the fraction (about one third) of recombinations within the bound atmosphere 
that go to the ground state and hence regenerate an ionizing photon. Hollenbach 
et al. ED presented analytic and approximate numerical solutions to the radiative 
transfer problem that imply ;?o(r) r“ 5 / 2 , and this scaling has been widely adopted 

in subsequent work. Important features of the Hollenbach et al. solution have been 
verified in more detailed radiation hydrodynamic simulations (including the <J> 1 2 
of the mass loss with the ionizing photon flux 112931 ), though the slope of no(r) at 
r > r g remains to be confirmed (indeed, a recent radiative transfer calculation is 
inconsistent with the canonical slope 032411 ). 
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8.1.2 Drivers of photoevaporation 

Photoevaporation driven by a photon flux of EUV radiation from the central star <f> 
is estimated to result in a mass loss rate mm. 



(253) 


provided that the disk is present and optically thick at all radii. The dominant uncer¬ 
tainty in applying this estimate to specific systems comes from lack of knowledge of 
<J>, which can be constrained but which remains hard to pin down precisely 15112651. 
For low mass stars (M < M 0 ) reasonable estimates imply that EUV photoevapo¬ 
ration rates are negligible when compared to the median accretion rate of T Tauri 
stars am but large enough to matter for disk dispersal if no stronger mass loss 
processes are operative. (EUV photon fluxes are of course vastly larger for massive 
stars, for which the theory was originally developed.) 

Low mass pre-main-sequence stars are strong emitters of FUV and X-ray radi¬ 
ation II1011274111321 1. The FUV luminosity has a base level that is set by chromo¬ 
spheric activity, on top of which there is a potentially much larger component from 
accretion cm The X-ray luminosity scales linearly with the bolometric luminos¬ 
ity, log(Lx/Lboi) = —3.6 Ii 274l . but with a large scatter. Qualitatively, these photons 
affect the disk in the same way as the EUV. The surfaces of the disk are heated, 
albeit to a somewhat lower temperature than the 10 4 K that is characteristic of HII 
regions, and where this heating results in unbound gas a wind ensues. Quantita¬ 
tively, the main difference is that X-ray and FUV heated layers are not separated 
from the cool underlying disk by any analog of a sharp ionization front, and this 
makes modeling of FUV and X-ray photoevaporation more difficult. State of the art 
calculations H128 .1 25611255 1. however, suggest that X-rays and FUV radiation drive 
mass loss rates that are substantially higher than the EUV prediction — with values 
of the order of 10 s M 0 yr 1 being possible — with X-rays likely dominating in 
the inner region. 

8.1.3 Disk evolution including photoevaporation 

Including photoevaporation in classical viscous models for disk evolution is partic¬ 
ularly simple, because thermal winds exert no torque on the disk (equation[63|. The 
rate and radial dependence of the mass loss, moreover, is primarily a function of the 
spectral energy distribution of the irradiation. This may depend upon the stellar ac¬ 
cretion rate (if the FUV luminosity is an important driver of photoevaporation) but 
it is not coupled at leading order to details of the disk structurj* 4 ] These properties 
mean that disks evolving under the joint action of viscosity and photoevaporation 


14 In more detail, however, the grain population within the disk will affect the absorption of high 
energy photons and hence the local mass loss rate fl29l - 
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Fig. 33 An illustrative calculation of disk evolution including photoevaporative mass loss. The 
model plotted is based on a disk with v °= r, and wind mass loss that scales as E w °= r outside 
5 AU. The disk displays the two time scale evolutionary behavior that is reasonably generic to 
internal photoevaporation models, with a long period of slow “viscous” evolution being followed 
by rapid inside-out dispersal. 


exhibit two distinct phases of evolution, an early phase in which viscosity dominates 
and a short subsequent phase in which the wind results in rapid disk dispersal f80l . 

Figure[33] based on the original calculations by Clarke et al. Il80l . shows how a 
disk evolves under the action of viscosity and EUV photoevaporation from the cen¬ 
tral star. In addition to the two time scale evolutionary behavior, the concentration 
of EUV mass loss toward the inner disk leads to a characteristic radial structure of 
disk dispersal. As the disk accretion rate drops, photoevaporation first dominates the 
evolution near to the innermost radius where mass loss is possible (in the illustrated 
example, this is taken to be 5 AU). A gap opens at this location, separating the inner 
disk (with a short viscous time scale) from the outer disk (where the viscous time 
scale is much longer). The inner disk then drains on to the star, stellar accretion 
ceases, and the disk is dispersed from the inside-out. The final dispersal is rapid, 
because the formation of a hole in the disk allows EUV (6) or X-ray 112571 radiation 
to directly illuminate the inner edge of the disk, which typically accelerates mass 
loss. The dust to gas ratio in whatever is left of the disk increases during the dis¬ 
persal phase 1327. 4l 1291. Although the details depend on the radial profile of the 
photoevaporative wind, broadly similar evolution is predicted in most models 11271 . 
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The observed lifetimes of disks are broadly consistent with theoretically estimates 
that are based on “best-guess” values of photoevaporative mass loss rates l22l . 

There is no universal form describing how the action of photoevaporation cuts off 
accretion in the inner disk. For simple models of EUV photoevaporation, however, 
it can be shown that the inner accretion rate M(t) is related to the accretion rate Mss 
that would be predicted by a self-similar model 112191 without mass loss via 12981 . 


M = 



(254) 


where to is the time at which accretion ceases. This formula is derived for a specific 
viscosity law (v«r) and photoevaporation model, but provides a qualitative idea of 
how the inner disk drains under more general circumstances. 

The inside-out character of photoevaporative dispersal applies only in the limit 
where radiation from the central star is dominant. In sufficiently rich stellar clusters, 
photoevaporation due to intense FUV radiation fields from other (massive) stars is 
more important. Unlike in the case of central star photoevaporation, for which the 
observational evidence is indirect 0 , photoevaporative flows driven by external UV 
fields can sometimes be seen directly, most spectacularly in the core of the Orion 
Nebula l42l . Adams et al. ni and Clarke 11771 have modeled the evolution of viscous 
disks under external photoevaporation, and shown that it results in destruction of 
the disk from the outside-in. For the fraction of stars that form in such clusters, 
this process evidently limits the time over which gas-rich disks would be present on 
scales comparable to the Solar System’s Kuiper Belt. 


8.2 MHD winds 


The theoretical and observational arguments for photoevaporation being an impor¬ 
tant component of disk evolution are strong, but it is not proven that other processes 
do not also contribute to disk dispersal. The obvious alternate candidate is MHD 
winds, which are likely to be present if mature disks retain a dynamically signifi¬ 
cant net magnetic field (see the discussion in §4.4.4 1. MHD disk winds have obvious 
qualitative differences from their photoevaporative cousins, 


• Their strength depends upon the disk’s net flux, rather than on the stellar radiation 
field. 

• They can be accelerated from arbitrarily small radii, where even EUV-ionized 
gas would be bound, with a velocity proportional to the Keplerian velocity at the 
magnetic field footpoint. 

• The local mass loss rate is (roughly) expected to scale with the disk surface den¬ 
sity, rather than being (approximately) a constant independent of the underlying 
column. 

• Predominantly neutral or molecular gas can, at least in principle, be accelerated 
(though it might subsequently be dissociated or ionized by stellar radiation). 
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If MHD disk winds are sufficiently strong, they can affect disk dispersal via some 
combination of mass and angular momentum loss. The resultant evolution can be 
quite similar to the photoevaporative case. In particular, if mass rather than angular 
momentum loss is dominant, Suzuki et al. BUI showed that MHD winds lead to 
the formation of a shallow surface density profile at small radii and, eventually, 
an expanding inner hole. In the opposite limit where angular momentum loss is 
strong (and mass loss negligible) Armitage et al. 03 suggested that dispersal could 
occur through the late onset of magnetic braking. Whether this is possible depends 
entirely on how the mass to flux ratio of the disk changes over time, and hence on 
the uncertain question of how net flux is transported and lost ((3.5.1 1 . 


It is highly likely that some level of net magnetic field threads disks whose sur¬ 
faces are heated to c s > vk by high energy radiation, leading to hybrid MHD / 
photoevaporative winds. There has been some work studying how MHD and 
line-driven winds interact f280l , but rather little is known of the observational 
or theoretical consequences of hybrid wind models in the protoplanetary disk 
context. 
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